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THEORY OF TRANSPORT IN LINEAR BIOLOGICAL SYSTEMS: 
Il. MULTIFLUX PROBLEMS 


JOHN L. STEPHENSON 
LABORATORY OF TECHNICAL DEVELOPMENT, NATIONAL HEART 
INSTITUTE, BETHESDA, MARYLAND 


If in a multiflux system the ith flux is given by the integral equation, 
t 

Y; = =, [ Ws; (¢— @)y; (@)da@ +m, (t), the corresponding equation in the 
0 


Laplace transforms is I; =2,;W;;I\;+M,;—the entire system having the 
matrix formulation, [1 -— W|I" =M. The general solution of this equation 
and its physical interpretation are discussed. Explicit solutions are given 
for the general mammillary and catenary systems and for some capillary 
exchange problems. The theory is applied to the integrated form of the 
fundamental continuity equation to give equations for total quantity of 
material in the various ‘‘compartments.’’ If the compartments are uni- 
formly mixed, the integral equation treatment is shown to be mathemati- 
cally equivalent to the usual differential equation formulation. 


1. Introduction 


In a previous paper (Stephenson, 1960, hereafter referred to as I), 
we discussed transport in a simple linear biological system, in 
which a single influx, y, (¢), and a single efflux, y, (¢), were re- 
lated by the integral equation, 


r= [ Pine a ado. (1) 
0 


In this paper we will treat multiflux systems which are bounded by 
several surfaces through which particles can pass in either direc- 
tion and which may be re-entrant. The analysis becomes more in- 
volved, but the underlying ideas remain the same. 
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2. Integral Equations for Multiflue System 


Before beginning the detailed and general mathematical discus- 
sion of our problem, a concrete illustration may be helpful. Let us 
suppose that a molecule of an isotopically labeled chemical com- 
pound (i.e., a ‘‘labeled particle’’) is introduced (‘‘injected’’) from 
the outside world into our system through a particular element of 
some surface. (This surface is a geometrical entity. It may or may 
not coincide with a physical surface such as a biological mem- 
brane.) Since the ultimate fate of this particle may depend on its 
vector velocity, we will suppose this lies in some infinitesimal 
range. This particle then wanders around in the system undergoing 
collisions with unlabeled particles. Collisions between labeled 
particles are assumed to be of negligible frequency, and the state 
of the system is assumed unaltered by the introduction of label. 
The densities in phase space of unlabeled particles are assumed to 
be stationary. In some of these collisions label may be transferred 
from the original molecule to another, or the particle may be con- 
verted into another species. Label, however, is conserved. Ulti- 
mately either the original labeled molecule or one which results 
from a collision process will pass through an element of some one 
of the surfaces with its vector velocity in some infinitesimal range. 
If we continue to observe the history of our label, such passages 
will occur repeatedly until the label passes out of the system to 
the outside world. The mathematical problem we will now attack 
is the description of the transport of label in such a system by a 
finite system of integral equations. 

We will assume that we have n surfaces of finite extent and m 
species of particles. Initially we will place no restrictions on the 
topological relations of the surfaces or on chemical conversions 
between species. Thus the surfaces may or may not be closed or 
parts of closed surfaces. If for the moment we consider only a 
single chemical species, a particle passing through any surface in 
either a positive or negative direction can newt pass through any 
other surface or the same surface in a positive or negative direc- 
tion. If we now consider all species and place no restrictions on 
chemical conversions (i.e., transfer of label from one species to 
another), a particle of a given species passing through one surface 
in a positive or negative direction can next contribute to the posi- 
tive or negative flux of any species through any surface, including 
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flux of the same species through the same surface in the same 
direction. 

In particular systems the probability of certain flux events oc- 
curing in sequence may be zero. Thus, if surfaces are mutually 
part of a closed surface, outgoing flux cannot contribute to itself 
but only to ingoing flux and vice versa. Likewise certain chemical 
conversions may not occur. When these physical restrictions arise 
they correspond to setting certain of the transfer functions equal 
to zero. 

If we classify our fluxes according to species, surface, and a 
positive or negative direction, we will then have a total of 2 mn 
fluxes, each of which can contribute to itself and all others. Our 
specific object is to relate these fluxes by a system of integral 
equations in which the contribution of the jth flux, y,(¢), to the zth 
flux, y, (¢), depends on a transfer function, w,,(¢- w). These equa- 
tions could of course simply be assumed, and we would have as in 
(7) below, 


y,(t) =, [ y, (0) w,;(t- 0)do + m,(t), 


0 


where m,(t) is a function accounting for the fraction of particles in 
y,(¢) not previously in some other flux, that is particles being in- 
troduced from some unspecified extraneous source. This approach, 
however, embodies some assumptions about microscopic events in 
the system which it seems worthwhile to examine in detail. In 
particular the subsequent history of a particle of a given species 
passing through a given surface will depend in general on the par- 
ticular infinitesimal element of surface which it passes through 
and on the exact direction and magnitude of its vector velocity. 
Clearly, to get a finite system of equations we must sum over finite 
surfaces and finite ranges of velocity. In doing this, .we will take 
as our basic assumption that for infinitesimal surface elements and 
infinitesimal ranges of velocity, flux is linearly related by trans- 
port functions such as the w,,(t - @) above, and will examine the 
problem of the four-fold summation over two velocity ranges and 
two surfaces to obtain the relation between two fluxes each of 
finite velocity range and through finite surface. As we shall see, 
two of the integrations are trivial and two are not. To define the 
transport functions relating infinitesimal velocity ranges and in- 
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finitesimal surfaces, we will suppose that a swarm of labeled par- 
ticles of species g, with velocity falling in the range A v, are “‘in- 
jected’”’ through surface element Ao,, belonging to finite surface S, 
at time w, the injection being instantaneous. Subsequent to this 
initial introduction of material into the system all labeled particles 
passing through any surface of interest will be captured, counted, 
identified as to species, and permanently removed from the system 
(by an appropriate demon). The flux of labeled particles of species 
k, passing through surface element Ao,, belonging to surface Si 
with velocity in the range Av,, expressed as the fraction of par- 
ticles initially injected will define our transport function. That is, 
we have 


Yijke ir "j> Ver %g1 ®t) Ao, Ao; Av, Av, AwAt = (2) 
Ly; (ris Ves o) Na, Av, Aallw;,; p, (Fists Ves v,»¢—w@)Aa;Av, At]; 


where Yijke (r;, r9 Vas V gs @, lb) Aa; Ao; Av, A v,AwAt is the number 
of particles of the kth species the velocities of which lie in the 
interval Av, which pass through the surface element Ag; with 
center coordinates r; in the time interval A ¢—these particles origi- 
nating from the ‘“‘injection’’ in our hypothetical experiment of the 
number of particles, Vie (r;,v¥,,@)Aa,Av, Ao, of the gth species 
with velocities lying in the range Av, through the surface element 
Ao, at time w. The second term on the right, Oi ke ay 
t-w)AtAo, Av,, is then defined by equation (2). It is the proba- 
bility that a particle of the gth species which passes through sur- 
face element Ago, with velocity in the range Av,, during the time 
interval Aw causes as the nevt flux event in the system the pas- 
sage through Ao, of a particle of the kth species, with vector veloc- 
ity lying in the range Av,, during the time interval Ag. 

As stated above, we assume that for infinitesimal surface ele- 
ments and velocity ranges the system is linear. Under this assump- 
tion, we can integrate (2) with respect to time to yield an integral 
equation. However, for any practical value in a system with sur- 
faces of finite area and with particle velocities distributed over 
finite ranges, equation (2) must be integrated not only over w, but 
also over rey Uys Vp, and v,- Note that integration over r; and rj 
corresponds to integration over S, and S;. Integration over r; and y, 
can be carried out with no difficulty (and the A t’s can be cancelled), 
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leading to an equation 


Yij,ke ("j9V p91 rt) Ao, Av, Aw = 5 

ly;, (r,,v,,@)Ao,Av, A ollw,; pe (tjr¥,st- a], an 
where yi; 4¢(Fjs¥g» ,¢)Aa;Av, Aw now refers to flux through a 
finite surface S, and a finite velocity range; likewise w,, ,, refers 
to transfer from infinitesimal surface Ag, to finite surface S, and 
from infinitesimal velocity range Av, to finite range of v,. The ex- 
perimental counterpart of this integration is that obviously after 
introduction of particles into the system through a particular sur- 
face element and within a particular infinitesimal velocity range, 
subsequently particles with velocities in any range and passing 
through a surface of any extent can be counted and summed. How- 
ever, integration of (3) over r, and v, is not so simple and in gen- 
eral will not lead to an expression of the form y,, (w) w;; 4, (t-w)Ao, 
but to one of the form ¢d(@, £-— w)Aw. In order to get a functional 
relation of the first form, either the transport function, w,;,, (rj, Veo 
t — w), must be independent of r, and v, over a finite range of each, 
or the distribution in flux of particles of the gth species passing 
through the jth surface must be separable with respect to time, i.e., 


Yig (r;5 ¥ 59 ®) = Vig (0) fj, (rj v,)> (4) 


In either of these cases we can carry out the integration. In doing 
this we will integrate over the entire range of variation of r,, but 
will separate the integration over velocities into separate terms for 
positive and negative flux.* A particle contributes to the positive 
flux or to the negative flux according to whether its velocity vector 
makes an angle greater or less than 7/2 with the positive normal to 
the surface element through which it passes. Carrying out these 
integrations and also the integration over previous time, we obtain 


*It should perhaps be additionally emphasized that the essence of the 
analysis in this section is the replacement of continuous distributions of 
flux in spatial and momentum coordinates with a finite number of discrete 
groups. It is obvious that the selection of these groups is to a large ex- 
tent arbitrary. For example, instead of dividing the velocity range into 
the two categories of positive and negative flux, it could equally well 
have been subdivided into multiple categories according to the magnitude 
and direction of the vector velocity. 
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the equation, 
t 
eitnee = | Yeje()Creiag(t- odo, (5) 
0 


where Y+i+;,k¢ (4) is the positive flux of particles of the kth species 
through the zth surface which is contributed by all previous posi- 
tive flux of particles of the gth species through the jth surface. To 
get the total positive flux of the kth species through the ith sur- 
face, we must sum over all surfaces (including the ith) and all 
species for both positive and negative flux; we must also add a 
term, m, ,,(¢), which will account for particles introduced from some 
extraneous source (i.e., particles that have not previously passed 
through any surface in the system); this we will call an ‘‘injec- 
tion’’ function. We then have for the total positive flux of the kth 
species through the zth surface 


¥+ ix (2) = DF Di Veiepke (¢) + oF 2, ¥+i-j,n¢(4) se m4 ;1,(t) + (6) 


where each of the 2nm fluxes on the right, some of which may be 
zero because of constraints, is given by a convolution integral of 
the type of (5). Before discussing the solution of this system 
of integral equations, we will introduce a simpler system of sub- 
scripts by reordering the fluxes. To do this we will order the 
fluxes, lexicographically, first according to surfaces and then ac- 
cording to species, identifying positive fluxes by odd number 
and negative fluxes by even numbers. According to this order- 
ing scheme, Y+j;g Will become y,, where A = 2(j -1)m+2(g - 1) +1, 
and y_,, = ¥,+,- From (6) we then obtain 


nO =%, [ w,(t- 0) y;(a)da+m,(t), (i =1,2..2nm), (7) 
0 


the subscripts on both fluxes and transfer functions referring to the 


new ordering. Taking the Laplace transforms of equations (7), we 
obtain 


DP, =2,W,0,+M,, (= 1, 2,..2nm). (8) 
These equations can be written in matrix form 


r=Wr+M, (9) 
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where I is a column matrix, the row transpose of which is 


RUNS eon Ions); (10) 
M is also a column matrix or vector, the transpose of which is 
Mie AM. =: «Mo, ..)s (11) 


and W is a square matrix, the entry of which in the zth row and jth 
column is W,,. In the following section solutions of (9) will be 
discussed. In this paper we adopt the point of view that the m, (¢) 
and w,,(¢) are known and seek solutions of the y,(¢) in terms of 
them. Frequently we are interested in the converse problem of de- 
termining the w,,(¢) from the y,(¢) and the m,(¢), but the treatment 
of this problem is postponed. 


3. Solution of Equations for Multiflux System 
A. General. Equation (9) can be written in the alternative form 
(d-W)r=M, (12) 


where | is the identity matrix. If we premultiply both sides of (12) 
by the inverse of (I — W) we obtain 


r=(1- WM. (13) 
Expanding (I — W)*, we obtain 
T= (Rew + Wee W.-M: (14) 


It will be noted that this expansion is in general convergent be- 
cause usually we have 


0 <W,,(p) £W,,(0) $1; (15) 
=, W,,(0) $1. (16) 


The terms on the right-hand side of (14) have a direct physical 
interpretation, the term W” M being the column matrix corresponding 
to those particles and only those particles making their (n + 1)th 
contribution to the various fluxes, which is to say to nth transitcon 
particles. This can be proved by induction: By definition the 
matrix M refers to those particles which are being introduced into 
the system and so making their initial contribution to some flux. 
The first transition particles in the ¢th flux are given by 2, W,;M;, 
which is the th term of the column matrix, WM. If we assume that 
the (n — 1)th transition particles correspond to the matrix, W"~ 1M, then 
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the nth transition particles in the ith flux are given by 
=W:; (W"~! M)., which is the th row of the column matrix, W (W"-1M) = 
W"M. Hence, the statement is proved for all n. 

The statement can also be demonstrated by direct construction: 
The component of the zth flux composed of (n + 1)th transition par- 
ticles, Yan 1(¢), is obtained by summing over all (n + 1)thtransitions 
which contribute to the 7th flux, that is, 


t w, 
y= E33, | my (¢~ 01) [ wih (iecorsYee 
y 0 0 (17) 
| . W, ;(@,-1 — ®,) m(o,)do,da,.. do,. 
0 


The Laplace transform of (17) is 
niles. 
U; Pe hee SLEW Wig oe Way My (18) 


which is the 7th term of W"M. 

B. Matrix partition. Frequently it is desirable to eliminate some 
of the fluxes from the system (8). Suppose the fluxes are divided 
into a set A, g in number, and a set B, A in number, g + h=2 nm, 
with corresponding division of injection functions. We then have 
the matrix equations, 


Cin WT, ewe Pe (19) 
Vp = Weal, + Weplp + Mp. ae 


It will be observed that W,, governs transitions within the set A, 
W,, transitions from set B to set A, etc. In general W,, and Wap 
being square will be non-singular and will have unique inverses. 
Solving (20) for", we obtain, 


Pr, =(I1- Wea) | (Wea I, + Mp); (21) 
(21) substituted into (19) yields, 
Py = Wag 4 + Wap (I - Wen) Wea rR 


; (22) 
Wap (l-Wep) *Mp+M,. 

Equation (22) can be solved for I", to give 

Ty, =(-W, - Wag (lpn) Weyl * (23) 


[Wop (l— Weep) 'M, + M1. 
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The various terms in ce have a physical interpretation: We Bet 
note that W, ,(I -W,,) 'W pa» Which is a square matrix of order g”, 
is the transition matrix for those particles which having contributed 
to a flux in set A, next make any number of transitions in set B, 
and then finally again contribute to a flux in set A. This can be 
seen from the definition of Maa» th ms expansion of (I - as in 


(14), and the definition of Wap? is 


Win (I — Won) Wa, = Wap <_ ee ee Was (Wee emer o2) 


The term W,,W,, enumerates the transitions 4 —> RB —> A; 
Wanton): oa; the’ transitions “A — B —> B..:. B= A= 
A — B” — A. From this interpretation of W,, (I — Wenz) 'Wa4, for 
which we will introduce the abbreviation Wiz.pa> the interpretation 
of (23) follows directly. 


We have: 


U — Wa, - Wapeal ‘UWap (I - Weg) ' Mp + My] = 
+ (W,, + Wap pa) t--- We + Win pa) t--! (25) 
[Wap (I - Wop) Mp + M4); 
Win (l -Wep) ‘Mp + M4] is the totality of particles making their 
first contribution to flux in A, and (W,, + Wp p,)" L.-J is the totality 
of particles making their nth contribution to flux in A. This can be 


seen from the expansion of (W,, + Wi, »,)” and is easily proved by 
induction if we consider the product, 


(Was + Wap, pa) (Waa + Wap,pa)” (Wap (I - Wop) Mp +My). 

We have developed these ideas in such detail both because of 
the insight they afford into the matrix formulation and also because 
sometimes in physical systems with relatively simple connectivity, 
a straightforward enumeration of fluxes according to transitions is 
the simplest approach to a problem. 

C. General mammillary system. Certain systems with particular 
connectivity are of special interest. One such is the general mam- 
millary system in which one flux can contribute both to itself by re- 
entrance and to all the other fluxes; the other fluxes can contribute 
only to the one flux, not to themselves or to each other. In this 


case our equations become: 


r,=2W,,0,+M,, ( assumes all values) (26) 


Wd 


Vr =Wl yeh, - ¢4D) 


| 
War) 
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The solution of the system (26) is, 
LW M.+M, 


pon lta 
ta ’ 
1- 20M, 
W, Cow M. + M,) 


1 eee eR Set gee Be pare 
‘ 1— > W_W, : 


f Siam Fes > | 


(Note that (26) and (27) are scalar and not matrix equations.) The 
analogue of this system for the compartmental type of analysis is 
that of a central compartment which communicates with a number of 
peripheral compartments which do not communicate with each other. 

As an example of the application of equations (26) and (27) we 
will suppose we can measure the flux of some material in the left 
ventricular outflow, y,(¢), and in the renal inflow, y,(¢), and that 
the material is introduced into the system by injection into the left 
ventricle. The above equations then become 


Dam Wy yb t+ Wig Va + My 


(27) 


and 
D, oo Wo r, 


In this formulation we have split the ventricular outflow into the 
three disjoint (pairwise non-overlapping) classes of particles ac- 
cording to prior history: (1) has been through some organ other 
than the kidney during the previous circulation, -given by Weta: 
(2) has been through the kidney during the previous circulation, 
W,.1,, (8) has just been injected, M,. It is physically obvious, as 
was shown mathematically in the previous section, that a particular 
partition of a given flux into classes of particles with different 
prior histories is arbitrary if the classes are both exhaustive and 
disjoint. 

D. General catenary system. In this system, of great practical 
importance, the flux through each surface is separated into a posi- 
tive and negative component. The positive flux through any sur- 
face is contributed to only by the negative flux through the same 
surface and the positive flux through the preceding surface. The 
negative flux is contributed to by the positive flux through the same 
surface and the negative flux through the succeeding surface. The 
positive flux through the first surface is introduced from the out- 
side. With no loss of generality we can assume that the system is 


TRANSPORT IN BIOLOGICAL SYSTEMS 123 


not generally re-entrant so that negative flux through the first sur- 
face and positive flux through the last are permanently removed 
from the system. We are thus led to the system of equations:* 


l,(+)=M,(4), 

66) = 8, G06) 4f,Or 0), 
Pre=F,®%F,4+R, OFC), 

TP, (+) ~ F,_, (+) 2) ae PR ; (-) lr; (-), (28) 
r;(-) =~ Fie O <r R,(4) 0,4) ’ 


a Wee ee Ole) Re Ce 6. 0. 6. 6 0) 6,6) 818 6, Ore Sree. 8 


i (+)= F 1 (+) ieee. +) ’ 
En (-) = 0, 


where I’, (+) and I’; (-—) are respectively the transforms of the posi- 
tive and negative flux through the ith surface; F;(+) is the trans- 
form of the transport function which relates positive flux through 
the ith surface to positive flux through the (7 + 1)th surface; F’, (-) 
is the transform of the function which relates negative flux through 
the ith surface to negative flux through the (2 - 1)th surface; F, (+) 
is the transform of the function which relates positive to negative 
flux through the ith surface; and Ff, (-) is the transform of the func- 
tion which relates negative to positive flux through the zth sur- 
face. It should be noted that the functions F (+) and FR; (-) de- 
termine the fraction of flux reflected before passage through any 
preceding or succeeding surfaces. The system of equations (28) 
can be solved by matrix methods, but a more useful solution is in 
terms of continued fractions. To develop this we introduce func- 
tions R,, (+) and F,(+)3 Rin (+) is the transform of the transport 
function which gives the total reflection of flux through the zth sur- 
face, that is, flux which is reflected after any number of passages 
through subsequent surfaces, and F,(+) is the transform of the 
function which gives the flux through the nth or final surface in 


*Note that although equations (28) are scalar, the algebra is developed 
so as to apply to matrix equations, provided of course that an inverse of a 


singular matrix does not occur. 
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terms of flux through the 7th, i.e.: 
Tr; (=) ee Ein (+) r; (+) ’ 


(29) 
Di (+) = F,,4)0,4). 
In particular, we have, 
T= FiG@Ts(), set 
Mi (-)=F,,(+)0,(), 
and also, 
0, (+) = Fa) 4), at 


r, (-) == R,, (+) | a (+) . 
The problem is to develop expressions for f,,(+) and F_,(+) in 
terms of F,,(+), Fy. (+), Ry (+), Fy (4), F,(-), and R,(-). Substitut- 
ing from (31) into (28) we obtain 


r, (-) oa R, (+) Yr, (+) uy F, (-) R,, (+) I, (+) ’ 
I’, (+) = Fy (+), (4) + RP, (-) Ra, (+) T, (+). 
The system (32) can be solved to give 
D, (+) = [1 - R, (-) Ry, (+)? Fi(+)T,(),; 
I, (-) = {Ry (+) + F, (-) F,, (+) [1 - Ry (-) Ro, (+))"? F, (+)} I, (+). 
From (33) and (31) we also have 
Ty, (+) = Fa(+) 1 - RB, (-) Rg, 4) F, (4) 1, 4). (34) 
Comparison of (33) and (34) with (30) and (81) gives: 
Fii(+) = Fy (+) [1 - Ra (-)R, (DIF, (4), 
Ry (+) =R, (+) + F, (-) R,, (+) [1 - Ry (-) Ra, (+17 : Fr (+), (35) 
: R, (+) te Fy (-) [R,,(+)* a Ry (-))- : F; (+) ° 


Equations (35) give recurrence relations for computing R,, and Ps; 
thus 


Fe, (+) = Ry (+) + F;(-) [R,, (+) * + R,(-)1"* F, (+), 
R., (+) = Re, (+) + Fi44(-) [Ri 4 1,n (+)°'- Ri, 1 (~)1" F,(4) . 


Substituted back into (35), the relations (36) give R,, (+) developed 
as a general continued fraction. Likewise we have 


Fi (+) = ores (+) U1 — Ria (-)Ryy in(t))? F,(+). (37) 


(32) 


(36) 
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If we introduce the abbreviation, 


Q,; = [1- Ryo, (-)Ritr,.@) FG), (38) 
we arrive at the final formula, 
I (+) = QO, ---e 0,0, M, (+)- (39) 


The above equations in addition to being valid for matrices also 
will of course hold in the limit when the transform parameter, p, 
equals zero. Then the various algebraic quantities represent sim- 
ply probabilities of reflection or of forward transmission integrated 
over all time. In certain cases the above formulae can be greatly 
simplified. For example, if 


R,(-) = Ra(-) =---- R(-) = +--+ 0) = 9, (40) 
then (39) becomes 
TD (+) = Fi-sG4) eee eee F, (+) F, (+) 4,@)- (41) 


If we consider the relations (35) for p = 0, i.e., 
[FH peo = Fai (+ 0) -[ f.(+5 t)dt, (42) 
0 


which gives the total probability that a particle of the entrant flux 
passes through the entire system without being reflected back 
through the surface of entry, and introduce the further simplifying 
assumptions (for p = 0) 


F (+; 0) = F,(-, 0), (43) 
R,(-, 0) =1-F,(-,0)=1-F,G, 9%), (44) 
R,, (+, 0) =1- Fa (+, 9); (45) 
we obtain 
1 t o (46) 


eee EE eeeennneens + —————-_ - 
Bag t+; 0) Fi(+,0)  FraaG, 0) 


Equation (46) has particular application to various gas transport 
problems. These applications have been discussed elsewhere 
(Stephenson, 1953, 1957, 1958). 

E. Capillary exchange. The above theory will now be applied to 
the problem of exchange between the vascular and the extravascu- 
lar fluid compartments. First we will consider a single capillary, 
with velocity of flow v, length L, and transit time ¢) (for particles 
which do not exchange with the extravascular space). We will 
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divide the capillary and its surrounding extravascular space into n 
series segments, each of length AL. We will suppose that for each 
capillary segment efflux is only into its own extra-capillary space 
and the next capillary segment, and that the efflux of each extra- 
vascular segment is only into the corresponding capillary segment. 
For the first segment we will have the equations 


Pies (47) 
T= WoT, + Wo3 Ts, (48) 
DT, = WoT, , ey 
Y= Wl a has (50) 


where I’, is the transform referring to influx into the capillary seg- 
ment, equal to the transform of the injection function, M,, for the 
first segment; I, refers to flux from the capillary into its extra- 
vascular space; I, to flux from the extravascular space back into 
the capillary, and I’, to efflux from the capillary segment (which is 
the influx for the next capillary segment). The corresponding matrix 
equation is, 


T={W+M, (51) 
where, 
rT = aia Be BA By t- (52) 
0 0 0 0 
W. 0 W. 
W = ee: Wy, ee : ; (53) 
Wi, 0 W,, 0 
MT = (M,, 0, 0, 0). (54) 


If we partition the fluxes into a set A, (T,, T,), and a set B, a0 
I',), the corresponding matrix partition is: 


(55) 


W. 0 0 W, 
Wi, = a} W “| 
BA ’ BR 

0 0 Ws. 0 


if we approximate (1 -W,,)~! by (I+ Wap), ie., neglect all but 
first-order transitions, equation (23) becomes on substituting equa- 
tions (55), 


(Py, Ty)? = (My, Wy My + Waa Waa Woy M,)". (56) 
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It will be noted that the essential result in (56), namely 
Dy = Way + Wy3 ¥ 5001) My, (57) 


can be obtained more directly in this special problem by suppress- 
ing the term W,,1T in (48). The expression (W,, +W,3W3;.Wo,) in 
(57) corresponds to F’, in equation (41). Hence for n identical seg- 
ments the transform of the efflux from the last segment is 


Pout = Was + Vas Wyo Wo1)” My - (58) 


For (58) to be of practical computational value it must be cast into 
analytical form: If the time required for a particle which does not 
exchange to go through the capillary segment is At, = AL/v =¢,/n, 
and if the fractional exchange coefficient per unit time from capil- 
lary to extravascular space is a, then a fraction of material aA7Z, 
will pass from the capillary into the extravascular space, and a 
fraction (1 - a@Aé,) will go through without exchanging. We have, 
therefore, 


w,,(t) =(1-aAt,) 5(¢- Az), (59) 
where 5(¢ — AZ.) is a delta function, 
W,,=(1-aAt,)e Pe, (60) 
W,,=aAt,, (61) 
Woo =W(p), (62) 
an arbitrary function of p for the present, 
Rf en PAto, (63) 


Strictly W,,, as given by equation (61), is for particles which enter 
the extravascular space immediately after entering the capillary 
segment; and W,,, as given by (63), is for particles which spend 
time AZ¢, in the capillary space after re-entry from the extravascu- 
lar space. However, regardless of the subdivision of AZ, into time 
before leaving and time after re-entry, the product W,,W,.Wo1 will 
be W(p)aAt, 27 PAts, Substituting equations (59) through (63) into 
(58), we obtain, 


Gin 4a 
ae fs ( a =| Seed hae W (7) “| zs 


x [:- sar aioe) ePto M.. 
n | 


(64) 
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As n — o, equation (64) assumes the limiting form, 
Tout = My exp {Lat W(p) - (a + pl ty}. (65) 


For the very special case in which the extravascular compartment 
for each segment is uniformly mixed with fractional exchange co- 
efficient & per unit time, then we have 


W=k/(p +k), (66) 


and, 
Pout = M, exp - {la + p- ak/(p +k) ty}. (67) 


For sudden injection of amount M, at time ¢ = 0, inversion of (67) 
gives the following expression for the efflux: 


¥ (2) cue = M, L8(¢ - ty) + e * (tt) (akt,)* (¢ ~ t,) 
I, (2 (akt,)? (t - t,)13, 


where /,[..] is a modified Bessel function.* This result can also 
be obtained by setting up a differential equation and solving it 
(Sangren and Sheppard, 1953). The differential equation treatment 
is only applicable to very special cases, however. The great power 
of the generalized formula (65) is apparent when we consider not a 
single capillary, but flow through a capillary bed with a distribu- 
tion f(¢,) of flow times for a material which does not exchange with 
the extracapillary space. If we consider flow of a substance which 
does exchange and assume that the capillaries are similar except 
for flow times (that is, the fractional exchange coefficient a and 
the function W(p) are the same), we have for the Laplace trans- 
form, H (p), of the transport function exchangeable material: 


H (p) = ip f(t_je ta tea ite dt, fe (69) 
0 


A 
2 


(68) 


If we compare (69) with the transform of the transport function for 
material which does not exchange, 


Pip)= [" He.)e™Ptedt, 
0 
we obtain the relation, 
H (p)= Fla+p—aW(p)]. (70) 


; *See E. T. Whittaker and G. N. Watson (1943, p. 372) for definition and 
discussion of modified Bessel functions, and A. Erdelyi et al., Bateman 
Manuscript Project, for inversion formulae. 
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This result has been stated previously and its significance dis- 
cussed (Stephenson, 1958), but without detailed derivation from the 
general theory. 


4. Equations for Total Quantity of Material 


A. General Solution. Frequently it is of interest to know the 
total quantity of material of a particular species which is present 
in a particular volume or ‘‘compartment’’ of a system being studied. 
Molecules of a particular species can arise either by transport in 
through the surface bounding the volume or by chemical reactions 
within it. Here, the basic equation is the integrated form of the 
fundamental continuity equations (see I): 


S; Vv 


i 


where S, is the surface bounding the 7th volume V,, J, is the net 
vector flux of the kth species out of the bounding surface, do is 
an oriented surface element, s, is the source function for the kth 
species, and q,, is the total amount of the kth species in V;. In 


order to compute the surface integral -{ J, -do we will suppose 
Si 

that the volume V, can be enclosed in a finite number of surfaces of 

finite extent, and that the positive and negative fluxes through these 

surfaces satisfy a system of integral equations, which can be 

solved by the methods developed in the preceding section. To com- 


pute the source integral / s, dV, we will suppose that for an in- 
Vi 

finitesimal volume element dV,, the source function for conversion 

of one species into another is given by 


8, (t,t) dV, = Ey ki yr (ts (rs AV, (72) 


where c,(r, ¢) is the concentration of the Jth species at position r 
and time ¢, and &; ,,(r, ¢) is the fractional rate of conversion of the 
Ith species to the &th at position r and time ¢. Integration over the 
volume V, gives: 


i S;, (r, t) dV, = Li a (es t)e,(r, t) dV; . (73) 
V; : 


i es 
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If the k, ,,(r, ¢) are constant throughout the volume V;, they can be 
removed from under the integral sign to give 


i 
= 2h: pi) Gi ()- 
Substitution of (74) into (71) gives the equation 


Vik (-,7z)- Vik (+, ¢) + 21k 41) Vil (¢) = dq;,/at, (75) 


where y,,(-, ¢) is the ¢otal flux entering V, and y;, (+, ¢) is the 
total flux leaving V, (the plus and minus signs are taken with refer- 
ence to the positive normal to the surface bounding V,). If the 
kK, ,,(¢) are constant with respect to time, the transform of equation 
(75) is 


8, (r, t) dV, = iki [ c,(r, t) aV,, 
Vi (74) 


i 


Vig) — Vin +) + 214i 41 Qin = PP, — Vi, O); (76) 
there being nm similar equations in the various species and 
volumes. 

If q;,(0) =0, the fluxes, y,,(-, ¢) and y,,(+, ¢), and the corre- 
sponding transforms can be determined from the subsidiary system 
of integral equations. The necessity for this restriction arises be- 
cause of the fact that in the solution of the system of integral equa- 
tions we assume that all the particles in the system have been in- 
troduced through some known surface at some known time, i.e., be- 
long to one of the ‘‘injection’’ functions, m,,(¢). If initially the 
system is filled with a known number of particles the spatial dis- 
tribution of which is unknown, this is equivalent to an entrant flux 
the time integral of which is known but the temporal distribution of 
which is not. If the spatial distribution is known initially, it might 
be possible to construct a corresponding ‘‘injection’’ function for 
prior time, but in general this would only be possible if a complete 
solution of the basic partial differential transport equation (see I) 


— div J, (r, ¢) + s(r, ¢t) = de, (r, t)/dt (77) 


were known, in which case the entire approach developed in this 
paper would be superfluous. 

It should be noted that it is on this question of initial spatial 
distribution (and on that of linearity) that the controversy regarding 
the integral equation description of metabolizing systems has hinged. 
In the initial formulation by H. Branson (1946, 1947), the time course 
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of metabolization of material already in the metabolic pool was as- 
sumed to be the same as for material entering the metabolic pool, 
As has been pointed out (Hearon, 1953; Wijsman, 1953), this will 
only be true if the material is metabolized by a simple first-order 
reaction out of a uniformly mixed pool; i.e., if its ‘‘age’’ relative to 
time of entry into the pool is irrelevant. 

Under the assumption that the system is initially empty and under 
the general restrictions with regard to linearity, the matrix equation 
corresponding to (76) is 


r(-)-T4)+KQ=71Q, (78) 
or 


Q=[p! -K) '(T(-) -TQ)]. (79) 
The column matrix, [[(—) — '(+)], can be computed from the injec- 
tion matrix, M, by the methods developed above: Denoting I (-) by 
[,, [(+) by Tg, M(-) by M,, and assuming M(+) = Mz = 0, we ob- 
tain from (21) and (23), 
Tp =(1- Wee) Weal ys (80) 
and 
Bee r= 
ee | Wes) | Weal ll — Nags Sap Al Wer) Weal My - (81) 
Substituting (81) into (79) we obtain finally, 
Q =[p!-K]"'[l ~(1W,,) * We] 
= Wa 4 — Wap Wes) pal 'My- (82) 


B. Equivalence of Differential and Integral Equation Formula- 
tion for Uniformly Mixed Compartments. The integral and differential 
equation formulation for uniformly mixed compartments must be 
formally equivalent if they are to give identical results. The demon- 
stration of this equivalence depends mostly on the proper definition 
of the entries in the W matrix and the utilization of some elementary 
matrix identities, particularly the identity, 


AB t= [BATI"'. (83) 


If we replace the surface integral in (71) by the approximation 
(see [), 


ah J, -do = Zhi, Gx (84) 
s 


i 
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we obtain the differential equation [cf. I, eq. (6)] 


Dis Tin + 21%, ar Ut = 19; / at, (85) 
with the transform 
mii Qin + Er hier Qu = PQin — Vix (0)- (86) 
The matrix equation corresponding to (86) is 
K°Q + KQ = 9!Q - q(0), (87) 


where the matrix K’ is constructed from the fractional turnover co- 
efficients Kies and K from the fractional conversion coefficients 
Ki gl in a manner which will be described in detail below. 

Initially, to avoid some complications in matrix algebra, we will 
consider the case in which there is a single chemical species. 
Here, K = 0, and (87) becomes 


(p! - K’)Q = q(0), (88) 
with the formal solution 
Q = (pl - K’)"* q(0). 
From (82) the corresponding integral equation formulation is 
Q = [pt U ~ (1 Wyn)? Weg] 
Ui Was — 4p (I - Wap)” Weal * M4. (89) 
Since no positive flux in the set A contributes directly to another 


positive flux, and no negative flux directly to another negative flux, 
we have the simplifying relations 


Waa = 9, Wap = 0; (90) 


and (89) becomes 
! 2 
Gs ae ~ Wealll - Wap Weyl * My. (91) 


The two column matrices or vectors M, and q(0) are identically 
equal;* hence, proving the equivalence of the two formulations re- 


*M, =q(0) may appear to contradict the above discussion on the 
necessity of the system being initially empty. This contradiction, how- 
ever, is only apparent because in that discussion q(0) refers to 
t=0-At=-—0, whereas in the differential equation q(0) refers to 
t=0+At= +0, where 0 designates the exact time at which introduction 
of material into the system is initiated. In this discussion of the equiva- 
lence of the two formulations (although not necessarily for the general 
integral equation treatment), it is assumed that the material is injected 
into each compartment as a bolus. 
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duces to proving the matrix identity 
= ! 
(PP Kt = <0 — Wy ~ Wag Waal”? (92) 


From (83), equation (92) is equivalent to 

pl—K’ = pl Wap Wp] l — Wy al? (93) 
To define the entries in the matrix W,,, which we will denote by 
Wpaij> We will suppose unit quantity of material is introduced into 


the 7th compartment at zero time and that all excreted material is 
permanently removed from the system. We then have 


q, (2) = exp (;, 2), (94) 


where k,, = - 2%,;%;;- Hence, for the efflux from the ith compart- 
ment we have 


vi (+, 8) = — ky exp (ke, 2); (95) 
from which it follows by definition 
Weait = —ky/(p — Kyi)- (96) 
Since only the 7th influx contributes to the 7th efflux, we also have 
Wpaiz = 9% res; (97) 


If we consider the efflux from the jth compartment, the fraction 
-k;;/k;; (which is a positive number) contributes to the influx of 
the ith compartment. Hence, we obtain the relations 


Wangan kaylkjp th, (98) 


and 
(99) 


W a Bii = 0. 


With W,, and W,, thus defined by equations (96), (97), (98), and 
(99), the identity (93) follows immediately by direct substitution 
and multiplication. * 

In discussing the general case of multiple species and multiple 
compartments, the first preliminary is to specify some ordering of 
the nm quantities, g,,, with a corresponding ordering of the fluxes 
y(+);, and y(-),,- This is necessary in order to actually write 
down the matrices K, K’, W,,, and Wp,. There are two natural 
orderings: One is to order first by compartments and then by species. 


*o[1—W le is a diagonal matrix with entries p —k;,; | — W,p Way is 
a matrix with diagonal entries 1 and off-diagonal entries - kj ae Ki;) 
Hence the product of the latter by the former is pIi-K. 
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Under this ordering we have 
QT = (O11) Q1a ++ Qims Garr +> Gams ++ Qnm): (100) 


The alternative ordering is first by species and then by compart- 
ments. Under this ordering we have 


=O 5 Oe let a Cy, 2 OL eats (101) 


The advantage of the first ordering is that it partitions K into n? 
matrices m x m each with only those along the principal diagonal 
different from zero. Thus, we have 


Eis = 
=>) K,, (102) 
i=1 


where K, is an m x m square matrix, bordered above by (7 — 1)m and 
below by (n — 2)m rows of zeros, and bordered on the left by (¢ — 1)m 
and on the right by (n —z)m columns of zeros. The entry in the 
kth row and /th column of the unbordered square matrix, K,, is Ke ere 
For matrices constructed in this way, we have the relations* 


K,K;i=0, ix}; (103) 


*Written out, we have for K,, K. and K: 


ki, 11 *i,12 te Meal ky 1m 
K; = ° . ; 
ki m1 ws ee. ig eee 
coe CO) ° 0 ee 
Bae Ae Ayepcts eee 
K,= Seer ‘ - Oh Ris Ps 
ee “imi ° *i em O.. 
ae UD . 0 Le Se oe 
K eee 0 


0 
O07 Ke 0 ef 


In general K, i 1 will exist. Therefore rea is K; ! bordered by zeros in the 


we onibs 
Same manner as Kye Note that K; “is $08 K; -1 , the latter not existing. 
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and =a 
Reka = E., (104) 


where E, is an m x m identity matrix bordered by zeros in the same 
manner as K,. Hence, we have 


2 ,E, =I, (105) 


where | is the nm x nm identity. From (103), (104), and (105) fol- 
lows the relation 


Se eae (106) 


With the second ordering K’ can be similarly constructed from m? 
square matrices, nxn each. Unfortunately for a given problem 
either one or the other orderings must be selected. Then either 
K or K’ will not have a simple construction. However, the transi- 
tion between the orderings can be made by introducing a permuta- 
tion matrix, which represents the transformation carrying Q into Q’, 
i.e., 


Q’=PQ, (107) 
and 
Q=P™'1Q’. (108) 
We then have 
K’ Q= P“1AQ’, (109) 


where A is the matrix constructed from the m matrices n x n each 
under the ordering Q’. Similar relations hold for A as the relations 
(102) through (106) which obtain for K. From (109) we obtain the 
following equations: 


KP AP: (110) 
K’-1=2PA1P-1, (111) 
A=PK’P™!, (112) 
A-t = P1K’~1P, (113) 


These relations give a convenient method of computing K’, since 
both A, P and P™! are easily found. However, frequently it is 
necessary only to indicate the transformation, as will be seen in 
the following discussion of the general case. 

The formal solution of (87) is 


Q=[p1=K - K’}"*4(0). (114) 
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For the corresponding solution in the integral equation formulation, 
we have, noting that as for a single species Wan = 0, W,, = 0: 


Q=([pl— Kis Lb W7 aw We ko (115) 
Factorization of (114) gives 
Q =[pl —- K]"*[I1 - K’ (pI — K)“']"1 q(+0). (116) 


Since as for a single species, M, = 4(+0), proving the equivalence 
of (115) and (116) reduces to demonstrating that 


Ul —K’(p1-K) "+ = (1 - We - Wye Welt (117) 


From (83), equation (117) is equivalent to 
LK (pt — K)0 SIE We Weg ih We ote oe (118) 


The proof of (118) depends essentially upon the definition of the 
matrices W,, and W,,, according to the fundamental scheme de- 
scribed in Section 2. In order to define Weai> We Consider the rela- 
tion between influx and efflux of the various species for the ith 
compartment, assuming that all particles are permanently removed 
from the system after their first contribution to efflux. If M,; is 
the injection matrix for the ith compartment, we have 


Q; = (pE, —- K, - Kj,) 'M,,, (119) 
where 


OF = (Qasr ee (120) 


and K‘,, is an m x m diagonal matrix with entry &;; , in the kth row 


n 
and kth column (Ais = — 24 Riek Ko:., being excretion entirely 
7=0 
out of the system), and E, is the mxm identity. If I',; is the 
matrix of the transforms of the effluxes of the various species from 
the 7th compartment, we have 


Tai = —K5; 9; (121) 
= — Ky, (pe; —Ki—R,) Mies (122) 
whence 


Weai = — Kp, (pE, -—K, -K4,)7}. (123) 
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If we partition W,, into m x m blocks, we have 


n n 
Woe ee i ys W455 (124) 

ii put 
where Weaij is the matrix in the 7th row and jth column of matrices 
properly bordered with rows and columns of m x m zero matrices. 
Since only the influx of the ith compartment contributes to the ef- 
flux, the off-diagonal blocks are zero, i.e., Weat = (0, 2#7; and 
Weaiz 1S Simply W,,,; appropriately bordered with zeros. It follows 

that 
Wea = — Kp (p1-K —Kj)°- (125) 
The matrix P= is constructed from the fundamental relations 


Paik = =Wasir Upy, > 7 At, (126) 


where WaBij, k =k;; ,/k;;,,- It is apparent that the ordering corre- 
sponding to (101) is the natural one for the construction of W,,. 
This is impossible since we have already selected the other for the 
construction of W,, and K. Therefore we utilize the transformation 
P defined by (107), and the relations (110), (111), (112), and (113) 
to give* 


W 


AR ed ih Bie Ap) (- Ap) *P 


= P~1[P(K’ — K*,)P~*][P(-K%,-')P"41P (127) 
= [K’ — K4][-K,]7. 


Introducing the abbreviation B = (pl —- K — K‘,) and utilizing (125) 
and (127) and (83), we obtain a ready proof of the equivalence of 
the two sides of (118), i.e. 


fa 4 = A yee 
H- 4. W,,)U—Wyal (=f -(K-K,) 8 ‘)U + KB] 
Sak Bt K? Bo 1]7! (128) 

=f—K‘(B + K;\"* 

=I—K’[p!-K]"!. 


*It should be noted that K’" = P-1A"P; also P =P ', since by defini- 
tion P represents the transformation carrying Qik into Qk; and so P* =/. 
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A model for storing information is constructed on the postulates that 
(i) traces are left within the model by threshold changes, (ii) these traces 
are widespread, (iii) the cells in the model have properties statistically 
determined, and (iv) the stream of information is itself the encoding and 
decoding device. 

The model properties are studied and are found to contain many of the 
basic phenomena associated with learning by living organisms. In par- 
ticular, this method of storing information exhibits a large storage capacity. 


1. Introduction. In the history of mathematical biophysics many 
authors have considered the phenomena associated with learning 
and have attempted to correlate them with the known properties of 
living organisms. As more information has accumulated concerning 
neurophysiological processes and as communication theory has de- 
veloped, it has become reasonably clear that certain basic principles 
are associated with learning as carried out by the brain. Any theory 
of learning must include such principles if it is to bear any rela- 
tion to reality. 

From the simplest conditioned reflex to the most complex be- 
havior of man associated with memory, there exists the common 
principle that past events have left an impression on the organism. 
In the case of the brain this principle has raised two related ques- 
tions: (i) does memory depend on a continuing activity or a static 
residue? and (ii) is the activity or the residue for each memory con- 
centrated in one area or distributed? 

Many nerve-net theories have depended upon memory or learning 
lying in continuing activity of neuron circuits. N. Rashevsky (1938, 
1945) developed this possibility and applied the idea of the rever- 
berating circuit with success to many aspects of learning theory. 
Again, W. S. McCulloch and W. Pitts (1943) showed that because of 
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the ‘‘all-or-none’’ character of nervous activity, neural events and 
the relations among them can be treated by means of propositional 
logic. A significant contribution by them (1947) has been their 
study of the perception of auditory and visual forms. 

There does exist, however, some evidence that memory, espe- 
cially of the long-term variety, depends upon physical traces being 
left in the nervous system involved. For example, experiments with 
hamsters show that the freezing of the animals or the withholding 
of oxygen until their brain rhythms cease does not damage their 
memories (Gerard, 1953). In their 1943 paper McCulloch and Pitts 
suggested that the phenomena of learning seemed to require the 
possibility of permanent alterations in the structure of nets and 
added that the simplest such alteration would be the formation of 
new synapses or equivalent local depressions of threshold. They 
suggested further that some axonal terminations, that at first were 
unable to excite the succeeding neuron, became synapses of the 
ordinary kind if the axonal terminations are excited at the same 
time as the firing of the neuron. Similar ideas have been put for- 
ward by many authors, including J. T. Culbertson (1950) and A. 
Shimbel (1950). Culbertson pointed out that threshold-variation 
models may be more economical in the sense that such networks 
perform with fewer neurons the same tasks as other model networks. 

With respect to the second question, it is probable that, what- 
ever the method of storing information in the brain, the traces re- 
lated to a particular memory are not localized in general but dis- 
tributed widely. This statement is of course subject to qualifi- 
cation. Direct evidence in favor of this view comes from K. §, 
Lashley’s work (1926, 1929). In so far as simple associations are 
concerned in the rat, at least one-half of the cerebral cortex can 
be destroyed without any appreciable effect, though damage is 
severe if a sensory organ or part of the nervous system used by 
incoming impulses is damaged, Lashley himself remarked that the 
memory traces might be localized individually without conflicting 
with the main facts, provided there were many traces and that they 
were scattered widely over the cerebral cortex (1929). With an ef- 
ficient filing system, the advantages of such a method of storing 
information are obvious. The main advantage is a safeguarding of 
memories against localized damage. A possible advantage is that 
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the method allows superposition of traces with a subsequent ‘‘en- 
larging’’ of storage space. A possible disadvantage is that the 
inflow of information will interfere with and ultimately destroy pre- 
vious traces. 

This view is bound up with the problem of the part the detailed 
neuroanatomy of the brain plays. In some early models of nerve 
nets, the prediction of future states depended on a complete knowl- 
edge of the state of every unit in the net at some instant of time. 
The shortcomings of such nets were summarized by A. Shimbel and 
A. Rapoport (1948) in the paper in which they developed their prob- 
abilistic theory of neural nets. In particular they doubted that 
the genetic constitution of a zygote carries a sufficient amount of 
information to predetermine the growth of a highly complex nervous 
system consisting of a tremendous number of minutely specific 
neural structures. W. R. Ashby (1952) also raised this question 
and concluded that the genes were too few to specify every neu- 
ronic interconnection. He suggested that the genes fix permanently 
certain function rules, but do not interfere with the function rules 
in their detailed application to particular situations. Thus the 
traces that make up a particular message may differ from nervous 
system to nervous system in detail, though the information in the 
message remains invariant. A. Shimbel (1950) has published a 
learning theory based on the lowering of thresholds under certain 
conditions in two ‘‘random’’ net models and has shown that con- 
ditioned responses can be exhibited by such models. 

Where the input to a nervous system is a constant stream of in- 
formation, the problem of the encoding and decoding of such in- 
formation arises. Granted that (i) physical traces are left by the 
lowering of thresholds, (ii) such memory traces are widely dis- 
tributed over the storage system, and (iii) the storage system is 
statistical in the Shimbel-Rapoport sense; what governs the dis- 
tribution of such traces through the system and recalls the memory 
later? We suggest that the information that entered during the time 
interval (¢ — 0) to ¢, 0 being a constant, acts upon the information 
entering at instant ¢, and that at some future time the recall of a 
piece of information which was inserted at instant ¢ is obtained by 
feeding in as an ‘‘unlocking key”’ the information fed in during the 
time interval (¢ — 6) to ¢. This principle, together with the above 
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three rules, forms the basis of the model discussed below. In this 
paper, only the simplest model consistent with these postulates is 
treated. It is hoped to treat more complicated models in future 
papers. 


2. Properties of the basic unit. It is supposed that information 
is accepted by the model on the binary scale through one channel and 
ejected by another. Also let time be ‘‘quantized’’ in pulse lengths 
so that the entry of a pulse length into or out of the model occupies 
one unit of time. A message is then taken to occupy a certain 
number of pulse lengths, some of which contain one pulse each. It 
is further supposed, following Shimbel and Rapoport (1948), that 
the refractory period of the model neuron con be omitted. This can 
be done if the refractory period is less than the synaptic delay so 
that a neuron may fire at two successive instants of our quantized 
time. 

Then Figure 1 shows the basic unit of the model with necessary 
fibres, It is to be emphasized that in this paper, the model used 
is defined primarily with respect to the four major postulates given 
in the introduction and not to be a faithful copy of neurons or neuron 
complexes found in nature. 


FIGURE 1 


In Figure 1; C is the initial lead-in, A and B are synapses, / is 
the pulse input, O is the pulse output, F is a feedback fibre, and 
L is a threshold limiter (hereafter called a limiter), 

A limiter has the property that it connects the thresholds of the 
two synapses A and B of a unit store in such a way that the lowering 
of the z-function of one (defined later) results in the raising of the 
z-function of the other by the same amount. 
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A unit store is the basic unit in the present model, which con- 
sists of a number of unit stores arranged in a manner described in 
Section 3. The action of the store is as follows: 

(a) Pulses arrive at A and B, 
(b) A and/or B fire or do not fire, 


(c) Limiter L changes A and B’s z-functions or does not change 
them. 


In a sequence of events (a) precedes (b) which precedes (c), 
Each process is now described in turn. 

(a) The activity at A and B is governed by the composition of 
the message and the geometry of the store. Let EB be equal in 
length to EA but let FB be different in length from FA and let these 
lengths be such that pulse lengths from F arrive synchronously 
with pulse lengths from E at A and B. Since some of these pulse 
lengths passing into the model have a pulse in them, there exists 
a set of ‘‘events’’ that can occur at A and B. For example, in one 
instant of time at A and B we may have pulse lengths with no pulse 
in them arriving from E while from F a pulse length containing a 
pulse arrives at A and an empty pulse length arrives at B. If the 
numbers 0, 1, 2 refer to the number of pulses arriving at A and B at 
a certain time, then the following eight ‘‘events’’ are possible at 
A and B: 


DrFRRNrROOCOrF y& 
OrRnwNOrRrR OF O&O & 


A B 
The example given above is the “‘event’’ 1 0. 
(b) Both A and B have thresholds for firing. A (or B) will always 


fire if two pulses are present at it (i.e., one pulse from £, the 
other from F’) and may fire if one pulse is present, depending upon 
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its threshold at that time. The threshold-function 7 therefore 
takes two values 1 and 2. If T is 1 the neuron will fire if one or 
two pulses arrive at it; if T is 2 the neuron fires only if two pulses 


reach it. 
A function 2, is now defined for both A and B which can take 


2z) + 1 discrete values -z), -z, + 1, Pg flo. esis , -2, -1, 0, 
Loway wees » 2) — 2, ) —~1, a. It is supposed that if 
~%y $@4 S$ +29, A fires for 2 pulses, and if 


~%y $4 < Xp, A fires for 1 pulse, where 


X, is the critical state where firing commences for the arrival of 
one pulse. Synapse B is governed by a similar function @ pe 

Thus the threshold T is a step function of a variable x which 
can assume 2z, + 1 discrete values, viz.: 


in or Tp 
Ca eae te 
Xn or Xp 
Xo Xo +Xo 
FIGURE 2 


When A (or B) fires, it transmits a pulse towards H which has a 
one-pulse threshold. 
(c) The limiter L links the z-functions of A and 8B so that 


C4 + 2p =0. 


Now let L operate, shifting x, and &p One place, for the following 
‘fevents’’ only, namely 


A B 


Event (i) 2 1 
Event (ii) 1 2. 


where as Ae the numbers 1 and 2 refer to the numbers of simul- 
taneous pilses impinging on A and B. 

Event (i): #4 is lowered one place, unless already in the ex- 
treme —@, position, while @, is raised one place, unless already 
in the extreme +2, position. 
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Event (ii): 2, is raised one place, unless already in the extreme 
+, position, while 2, is lowered one place unless already in the 
extreme ~—2, position. 

This shift is termed a limiter shift and it must be emphasized 
that it occurs only when (i) E has sent pulses to A and B, and (ii) 
A or B but not both has a pulse from F. 

Thus facilitation is achieved since, if either A or B is fired 
often, its threshold tends to be lowered until it fires for 1 pulse 
only. But because of the relation 


C4 ™—%p> 


lowering the threshold of A may lead to increased threshold of B 
(depending upon the position of X,). C, the initial lead-in, is nec- 
essary in the unit store by itself, since in general it is possible 
that neither A nor B has a threshold low enough for one pulse to 
fire it. In the complete store described in the next section, C is 
not necessary and is omitted. 


3. Design of the information store. This store consists of a set 
of (n — 1) unit stores arranged as in Figure 3 where the symbol 


denotes a unit store. 


FIGURE 3 
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The neurons D,, D,,......, D,_, have synapses of unit threshold 
such that the arrival of one pulse at them causes them to fire, 
transmitting one pulse. The neurons H,, H,,......, H,_, are sim- 
ilar in character. 

G is a cell of fixed threshold such that it will fire if and only if 
a sufficient number of pulses reach it during one unit of time. The 
threshold is fixed such that this number is above the ‘‘noise 
level,’’ the average number of pulses reaching G and due to the 
random firing of the A and B synapses in the unit stores E,, 
E,, ..., ©, _) (see Sec. 4). 

Fis Eageaens F’_, are the feedbacks from G to the unit stores, 
each F; going to A; and 8; such that FA, is not equal to FB. 
Though drawn simply in the diagram they must include internun- 
cials acting as delay lines such that when a message is passing 
through the information store the feedback action is progressively 
delayed in going from E, to E,_, in the manner described below. 

Consider a stream of information entering the store. This stream 
consists of pulse lengths, some of which contain one pulse each. 
Let the stream enter the store at / at a given rate of so many pulse 
lengths per second. Label any consecutive set of (n +1) pulse 
lengths 1, 2, 3,......, n-1, n, n+ 1, in the order in which they 
will enter the store. The feedbacks F, to F_, have lengths M 
i.e., numbers of intermediate neurons M so chosen that when the 
(n + 1)th pulse length, duplicated by all the D; neurons (Fig. 3), 
has reached the A and B synapses of all unit stores by way of the 
D; neurons, the preceding n pulse lengths are distributed at the A 
and B synapses (having arrived by the feedback paths F, to F_,) 
as follows: 


In unit store £',, pulse lengths 1 and 2 have reached A, and B, 
respectively, 


In unit store F,, pulse lengths 2 and 3 have reached A, and B, 
respectively, 

In unit store E,, pulse lengths 3 and 4 have reached A, and B, 
respectively, and so on until in unit store Ey pulse lengths 
(n — 1) and n have reached A,_, and B,_, respectively. 

Thus when a message is entering the information store, each n 
pulse lengths of the message may be said to ‘‘act upon’? the 
(n +1)th pulse length and govern its passage through the unit 
stores according to the following argument. 

The arrival of a pulse length from the D; neurons at the A and B 
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synapses of the unit stores will in general precede the firing of a 
spectrum of pulses S by the unit stores via the H, neurons to the 
““collector’’ neuron G. It is seen later in the paper that at G, the 
pulse length that originated this activity is reconstituted and sent 
to O, i.e., a pulse arriving at / and then being sent simultaneously 
by the D, neurons to the A and B synapses is remade and sent by 
the G neuron to 0. The spectrum of pulses S can then be said to 


be the passage of this pulse through the unit stores between D; 
and H; But § is also a function of the particular spectrum of 
pulses of the preceding n pulse lengths of the message and there- 
fore each n pulse lengths of the message ‘‘act upon’’ the (n + 1)th 
pulse length and govern its passage through the unit stores. 

In the analysis that follows, a typical message is considered to 
be a stream of pulse lengths, some of which contain a pulse each. 
In practice, a language such as English is highly redundant so that 
its information content is less than a synthetic language made up 
of the letters of the English alphabet used with equal frequency. 
We assume here that there is no redundancy in the messages fed to 
the store, the probability of a pulse length being filled (i.e., con- 
taining one pulse) being } so that information is being presented 
to the store at a maximum rate on the binary scale. 

Again, since the probability of a pulse length containing a pulse 


is 3, the probability of obtaining in a set of n pulse lengths (n — 7) 


pulses is 
ae al r 
n@(5)" (5) <2". 


If the set is taken large enough, i.e., if n is large, the probability 
of any particular deviation from equality in numbers of empty pulse 
lengths and filled pulse lengths (i.e., containing a pulse) in a set 
of n pulse lengths falls off very rapidly with deviation. We take n 
to be so large that appreciable deviations from equality in numbers 
of pulses and pulse gaps do not occur so that we can consider 
average conditions in the store. Effects dependent upon devia- 
tions from average are considered later. 


4, The noise level in the store. Let a message be fed in to the 
store and let it end abruptly at a time such that the last n pulse 
lengths of the message are fed back to the unit stores and have ar- 
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rived there, but there is no pulse being fed in by the D. cells. 
Also choose this time such that the particular set of n pulse 
lengths has not previously been presented to the store when fol- 
lowed by a pulse in the (n + 1)th pulse length. 

These conditions are necessary for the following reasons. 

Even in the absence of an input to the A and B Synapses from 
the D; neurons, any feedback spectrum of pulses will trigger-off 
palsies from those A and B synapses that have at that moment one- 
pulse thresholds. It is therefore essential to consider what will 
happen, for example, at the end of a message, i.e., when its last n 
pulse lengths are fed back by the F, toF,_, paths. 

In addition, in order to see what effect any feedback set of n 
pulse lengths will have by itself, we must make the proviso that 
the set we choose in this section of the analysis has not pre- 
viously been a feedback set which operated when a pulse had 
reached the A and B synapses via the D, neurons. For if it had 
been, a ‘‘memory”’ of this feedback set would have been left in the 
thresholds’ distribution, as is shown in Sections 5 and 6, thus 
making this feedback set invalid as a test for ‘‘noise.’? 

Now each threshold limiter in a unit store can exist in one of 
(2¢, + 1) states, namely 


A B 
—2, +25 
—% +1 +% —1 
Sao +2 -—2 
: 22, + 1 states. 
+2) -1 —%,+1 
+2, —fp 


Let the (n ~1) limiters be Scattered uniformly throughout the 


threshold 2-function range so that there are 


1 
a limiters, on the 
average, in each state. With n being large, there will be from the 
meee feedback, 


“LA ite 
(ye cera limiters receiving 1 pulse at A, 0 pulse at B, 
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(ii) ~ eT ee receiving 0 pulse at A, 1 pulse at B, 
imiters receiving 1 pulse at A, 1 pulse at B, 

ee Set ee i. 

(iv) limiters receiving 0 pulse at A, 0 pulse at B. 


Then since there is no pulse being fed in by the D. neurons in 
this section of the analysis, the number of firings will be, from 


je. 


: -(z,) + X, + 1), from the A’s, 


Pd 
_—s 
hs 

— 


; -(%) + X,_ + 1), from the B’s, 


-~—-_-~ 
_. 
—_s 
—s 

— 


¥ (2) + X_ +1), of which 


—__——-. X,+1 from the A’ d 
SSG (2. + X, + 1) are from the A’s an 


He] Ae] ee | DOR] eR] 
3 
! 


» (2% + X, + 1) are from the B’s, 


-~o~ 
nw, 
< 

_— 
oS 


%yo+X_,+1 


It may be noted in passing that - 
z 


is the probability of 
o + 
finding a value of z such that 
—%, <#%SXo, 

also that a unit store (see (iii) above) may in certain cases send 
two pulses to its H, neuron, if the feedback to it consists of two 
pulses and both A and B are in the one-pulse threshold region. 
This last property depends upon the value of X, and may exist for 
the model X, 2 0. 

Referring now to Figure 3, remembering that the H; neurons have 
one-pulse thresholds, and considering the two cases X, <0 and 
X, 20 separately, it is seen that the number of pulses fired to- 


wards G from the H, neurons will be 


ai n—- 


1 
Eig mae sige Xp 20, 


(i) 


(iy 7 ae + (2% + Xo IE <0, X, 2.0, 
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Rot iad One 

OE Baer 1k On NOE eaten 
1 n-1 
UE EE 7 1), Kea 
4° Oe ey (Soated) ceo 


(iv) 0, Xo = 0, sp 20, 


giving, in total, 


n-1 

Se. a (#0 + Xo +1), for X, <0, 
0 

n-1 Ae 

Sey si’ (ot ght g) or Xoo 


These numbers may be taken as noise levels since at any in- 
stant, as we have seen, any set of n pulse lengths used as a feed- 
back will trigger-off a number of pulses from the A and B synapses 
even if no input pulse is entering via the D, neurons at that in- 
stant, and if this particular set of n pulse lengths has not pre- 
viously been used as a feedback when an input pulse entered via 
the D; neurons, the limiters will not have changed thus no record of 
this set will exist in the store (see Sec. 5). The output of pulses 
sent by this set to the neuron @ is consequently not significant 
and is referred to as “‘noise.’? The average noise level is taken to 
be the above-defined noise level when the set of n pulse lengths 


n 
contains ; Pulses and °) pulse gaps. 


The threshold of G is taken to be hk, where 


———— + (4%) +X, + 1) <h, Ay <0. 


Say 8 <h Xn 0 
Deviations from average are considered in Section 8. 


5. The storing of information. Let a message enter the store 
and let us consider any length of it consisting of (n + 1) pulse 
lengths. There are two cases to be considered: 

(a) where the (n + 1)th pulse length contains no pulse, 
(b) where the (n + 1)th pulse length contains a pulse. 

Case (a): Since no limiters change position, no change is made 

in the set of limiter positions. Also, 
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n-1 


22, 1 


-(%5 + 45: + 1), for X5 <0, 
or 
n-1 


KG 3 for X 
Qe, +1 Polo a adie toe 0 20, 


pulses, the noise level, are fired to the neuron G which-does not 


fire. When the first n pulse lengths of this section of the message 
are reinserted in the store, 


n-1 
. X 
on. + 1 (a, +A, + 1), for X, <0, 
or 
n-1 Ayiics 
ilecciies Zt +7) for X, 20, 


pulses, as before, are triggered-off, giving no output past G. 
Case (6): On the average, there will be from the message feed- 
back, i to the first n pulse lengths, 


(i) “ eyeiiters receiving 1 pulse at A, 0 pulse at B, 
(ii) * eSiaiters receiving 0 pulse at A, 1 pulse at B, 
(iii) “ Myinitors receiving 1 pulse at A, 1 pulse at B, 
(iv) 2 I limiters receiving 0 pulse at A, 0 pulse at B. 


In Re all limiters receive at both A and B 1 pulse from the 
input (the (n + 1)th pulse). 


=a - ' 
Thus is limiters, those of cases (i) and (ii), change their 


1 an- 
state by one, except the —-+ 5 which are already in the ex- 


2 2, + 
treme positions -2,, +#,) and +a), — % and cannot move. 


1 
ticular, —- 
In particu aS Smlegit 


a < X, to within it iar a = X,) so that the number of limiters in this 
set Aiden by the pulses of the message that have one of their A 
n-1 


move from just outside the range -2%) < 


1 
or B synapses firing for 1 pulse is augmented by —. 
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It is seen that at this time, the number of pulses sent to the G 
neuron arises as follows: 
Case X, 2 0: From 


a1. m~1 2. +X +1 
z A synapses fire, also oki Eekabs Bale 


(i) B synapses fire, 


2¢,+1 


n-1 : 
resulting in pulses being sent from their H, neurons to G. 


- p Gere eb m—1 
Lod lon Soe A synapses fire, also 
Ve es! 


B synapses 


(11) 


n-1 : 
fire, resulting in pulses being sent from their H, neurons to G. 


n-1 n-1 


A synapses fire, also B synapses fire, resulting in 


(iii) 


n-1 


pulses being sent from their H, neurons to G. 


o— eX ed ij] se og: Neer 
-2__0 __ 4 synapses fire, also af EME be Soe: 


igre aprhanes 


nN — 
Synapses fire, resulting in pulses being sent from their H; 


neurons to G, 
Therefore, for the case X9 2.0, a total of (n-1) pulses reach 
the collector neuron G. 


nm—1 
Case X, < 0: From (i), (ii), (iii) as above, pulses reach G, 


From (iv) however, since X, <0, only 


n#=1T Slee Xe T 
nae Eo #0 *)) pulses 


22, +1 
reach G. 
Therefore, for the case X, <0, a total of 


m-1 8e,+2X +5 n-1 Xx A 
0 ee OF ot a ee 
4 Qe, +1 Qe, +1 ( ae eer 
pulses reach the collector neuron @. 
These numbers are taken to be sufficient to fire the cell G so 
that a message being presented to the store passes right through it. 
When the first n pulse lengths of the message are now reinserted 
in the store, the same feedback of pulses and pulse spaces is pre- 
sented to the array of limiters changed on the previous occasion. 
It must be considered now whether the pulse that existed in the 
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(n + 1)th pulse length is remembered or not; in other words, has the 
store the ability to reconstitute that particular pulse by the change 
it made in the limiters’ positions, to feedback this pulse with the 
previous (n~—1) pulse lengths, and continue to recall the rest of 
the message? 


6. The information store’s memory. Once again there are two 
cases to be considered: 


(a) No further information has been stored between the presentation 
of the message and the re-presentation of the message’s first n 
pulse lengths. For example, the ‘‘organism’’ has slept. 


(b) Other information, that does not include the particular set of 
(n +1) pulse lengths under discussion, has been fed in between 
the presentation of the message and the re-presentation of the first 
n pulse lengths of it. 

Case (a): The positions of the limiters are as they were left 
on the first presentation of the message. Now let the first n 
pulse lengths be fed in to give the same feedback of pulses and 
pulse-spaces. The first presentation of the message shifted 
Lt —~ 1 
2 Qe +1 
easily shown by the above arguments that the response of pulses 
sent to G is now 


n-1 3 
Se te fe oe Ce XK <0 
22, +1 (2, >): Br 


limiters into the one-pulse threshold range, and it is 


a yas 
“s a aa Xs 0: 
Qe, +1 2 


These numbers are greater than their respective noise levels and 
if h is chosen such that 


ast 3 
2 | +X +5) > A> 
Qe, +1 0 ) 


or 


- -1 X 3 
a RE +st)oar 2. eta X, 29, 
Qe, +1 ce ate 2¢, +1 e 


then neuron @ fires and the store ‘‘remembers”’ the ( + 1)th pulse. 
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It is obvious from the above arguments that if the message had 
been presented twice instead of once (including the (n + 1)th pulse 
both times), an additional number of synapses would have been put 
into the one-pulse threshold, showing that the more often a piece 
of information is presented, the more firmly it will be fixed. 

It is to be noted also that although a recall greater than noise 
level is achieved by the message being considered, any other set 
of n pulse lengths will, being different from the set under con- 
sideration, simply give on average a noise-level response since 
its pulses are as equally likely to select A, B synapses the 
thresholds of which have been raised above the one-pulse threshold 
as to select A, B synapses the thresholds of which have been 
lowered below the one-pulse threshold. It is only, in fact, mes- 
sages very similar to the set under discussion that have a chance 
of giving a response greater than noise level. In a sense, there- 
fore, the store ‘‘looks’’ for situations in its memory store similar 
to the situation with which it is faced. 

Case (b): Let time ¢ be measured from ¢ = 0, the time of presenta- 
tion of the original message. In the store there are (n — 1) limiters 
and on the average a pulse enters the input every two pulse lengths 
while, on the average, the number of limiters having feedbacks of 
either 1 pulse at A and 0 pulse at B or 0 pulse at A and 1 pulse 


. n - 
at B is Thus the probability of a particular limiter changing 


1 
at any instant is 7 This is not strictly so (see Sec. 8) but suf- 


fices at this stage. 

We choose the most pessimistic case, namely that after the 
original message of (n +1) pulse lengths has been stored once 
only, a message of length m pulse lengths, where m is much 
greater than n, no (n +1) consecutive pulse lengths of the new 
message being identical with the original message, is fed in con- 
tinuously. Now time ¢ is measured in pulse lengths so that 8¢ 
equals 1 pulse length. We require to know how long a message is 
necessary (i.e., how long m must be) before the original message 
is ‘blurred’? beyond recall so that its output of pulses to G when 
the first n pulse lengths of it are fed in falls below h. 

When the original message was stored, its first n pulse lengths 


limiters 


were fed back to the limiters so that on the average 
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n 
went one stage toward —2, for 4, 


limiters went one stage to- 


n-1 


ward —2#, for B, and limiters stayed put. 


The change in this specific arrangement of thresholds with re- 
spect to the test message only as the new message of indefinite 
length is fed in is what must now be studied. 

Before storing the test message the distribution of limiters is 
as shown in Figure 4 where the ordinate y is the number of limiters 
in each state. 


Y 

A+x, x oO =x. =xoA 

B-x, —X, O- X; +x,8 
FIGURE 4 

ae on Sere 
It is seen that there are ————— limiters in each state, 
w+ 
0 
A B A B 
to 
-2, +2, +2, -2); 


X_, being the one-pulse threshold. 
After storing the message of (n +1) pulse lengths, the distribu- 
tion is as shown in Figure 5. 


QGQAAAGGGWW 


WB 


A+Xo Xo 0 =Xo XA 
B-x, 2S 0 Xo +%_B 
FIGURE 5 
1 a”-l1 ; ay 
The unshaded set, of depth — -————,, consists of those limit- 
2 ee ee. 


ers unaffected by the message. 
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n * . . 
The cross-hatched set consists of the limiters sent one 


stage toward —2) for A and therefore has a distribution 


y = 0, for #=+2, 
Lene with respect 
= — +-—___ ~1>¢>- 
Yy A renee for @, -— bees a, +4 eS: 
1 ; z-axis of A, 
nr— 
SS 5 for z= -2 
so 2¢,+1° 


nN _ 
The stippled set consists of the 


limiters sent one stage 


toward —2, for B and therefore has a distribution 


y = 0, for e@=+2, 
1 a-1 : 
yu: ‘ for @,-l>e>-2,+1 with respect 
4 2¢, +1 to the 
5 ee t-axis of B. 
ae : for z=-2 
OO erat 0 


These three sets can therefore be combined conveniently as 
follows, but only with respect to the test message: 


FIGURE 6 


Here y, is the number of threshold e-functions in state z at time 


¢. Immediately after the test message has been stored, that is, at 
t=0, 


3 ana-1l ; 
‘sa ea ye a 
n-1 


’ for “e,+lse<¢a,-1 
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y, =— +———-, for t=+2). 


As ¢ tends to infinity, this profile is smoothed down to a uniform 
n-1 
2a 
of indefinite length is being fed in, change as follows at any in- 
stant of possible change (a pulse being fed in at J): one quarter of 
them move one step toward —2_; one quarter of them move one step 
toward +2,; one half of them remain as they are, except those in 
the states -2, and +2). Because they are in the extreme positions 
at the barriers, three quarters of them remain in this state while 
one quarter of them move into state -¢, +1 and state +z, -—1 re- 
spectively. Thus the distribution in Figure 6 is gradually smoothed 
out. 

Since the instants of change are on the average every two pulse 
lengths, though subject to Section 8’s refinement, the equations 
governing the number y at any time ¢ after the message has been 
stored are set up as shown below. 

Let Sy, be the change in the number of thresholds in state # in 
a time S¢. 

Then on the average, 


1 1 1 
by. - (Fy gr ) oe 
% 8 x+1 4°% 8 x—1 


dy, 1 1 1 
are gixti ~ ee ~ Buen ’ 


giving the set of equations, 


dy..A 1 1 

Ge atstt ghey FTE 
dy 1 1 1 
Cina: meee oe el 


depth 7 since the threshold z-functions, while the message 


or 


e=-2, (1) 


Lp EN tee aes 
Je = Hy — Ty +=y 


gre et es ee eta 


7. The solution of the limiters’ equations of motion. Before giving 
the solution to set (1) let us consider the relation of #, to X). 
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From Section 6 we had 


n-1 3 n-1 
(2, +x +5) > eg oth +d A, 


v-1 X 5 n-1 X 3 
& (2 tne) > as (« +24), AP ows 
Sc 1 tak Es oe 


If z, is large with X, positive and much less than z,, then / is 
confined to a very narrow region, that is, the signal-to-noise ratio 
is almost unity. If X,= ~2,, then 

3°) pat n-1 
ae > h>—__—__, 
2 2¢,+1 i ae 


and the signal-to-noise ratio is 1.5. The investigation will be con- 
tinued for this particular model. Certainly this ratio would be im- 
proved by running the message through more than once, but at 
present this feature will be ignored. 

Then, letting ¢ = 81, the set of equations (1) may be rewritten 


This set of linear first-order equations is of the form 


ay, oe ; 
pee hoe G.-Y, t=1,2,....., m, 


Hon 


where m = 2a » + 1 and the @;, are the appropriate constants), 


A complete discussion of this set of equations is given by E. 
Trucco (1960). 


The general solution is 


m-1 
k a 
Yg(t) = )° 18, cos lex -1 =| °é Ze O=1,2,...,m, 


k=0 
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kn : kn 
Ba (1 ~ cos 7) = 4 sin? (#), foek =O, 102) 2o mre ls 
2m 
The quantities B, are the m required arbitrary constants, to be de- 
termined from the initial conditions. These conditions at T = 0 
were 


¥_.. = So: 
2 
Y_ = 809 forsee ie 2, — 1 
1 
Y,. age 
where 
3 nm—-1 
Gee : 
2 22, +1 


Trucco (1960) shows that under these conditions the arbitrary 
constants B,, where & is even, vanish while the constants for odd 
k are given by 


Thus the solution finally becomes 


m= 1 
2 4 Bs 


Ihe! 
where 


A} 3m 2m 2m 


Qi-—1 
and the apie 4 sin? aoe 5 


Ba Chains a se: Pw 


It is seen that: 
(i) no y, tends to infinity; 
(ii) a constant term appears as expected (for & = 0), each y, tend- 
ing to it as T tends to infinity; 
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1 
(iii) in the exponential of smallest root B,, B, tends to Ee 


2 
7 


(Qe, +1)? 
states tends to infinity; 


as @, tends to infinity, that is, as the number of limiter 


(iv) in the exponential of largest root B ,-1» this root tends to the 
2 


value 4 as az, tends to infinity; 


(v) the constant C,, becomes smaller as #, increases, behaving 
approximately as (2¢, +1)’. 

In general, therefore, as T increases, the exponentials of small- 
est B, ~ (22), +1) 7, rapidly become the dominant ones so that 
essentially ee the number of pulses fired toward the cell G when 
the first n pulse lengths of the test message is reinserted, decays 

ima 


(2x9 + 1)" : 
finally as the exponential e » though more rapidly at first. 


As an approximation, therefore, remembering that ¢ = ah, 


: mates 
ie (2x5 + 1)° 

| (24.0, ; ) 
Thus, it appears reasunable to take U, given by 


1.23 U 
(Qe, +1)? 


= Le, Ora tt 


aS a measure of the ‘‘life’? of a memory in the store when the mes- 
Sage giving rise to the memory has been stored once only. 

Before proceeding further it is useful to summarize the argu- 
ments leading to the analysis of Sections 4,5, 6, and 7. It is seen 
that an important property of the store is that the sequence of 
events after n time units depends on which particular set of n 
pulse lengths is selected for analysis. 

It is seen that any set P of n pulse lengths arriving by feedback 
at the A and B synapses ata given instant is in one of two classes: 


(1) at that instant no pulse has reached the A and B Synapses from 
the input via the D; neurons, 


(2) at that instant a pulse has reached all the A and B synapses 
from the input via the D, neurons. 
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Class (1): The set P, by itself triggers-off a spectrum of pulses 
to the H, neurons and hence to the collector neuron G, and this 
spectrum is characteristic of that particular set P,- This class 
(1) is now divided into two subclasses (a) and (b), where in (a) the 
set P, has never previously been inserted in the store, and (b) the 
set P, was at least once previously part of a message inserted in 
the store. 

Subclass (1a): This set P,, is not significant to the store, i.e., 
it is not one of its ‘“‘memories.’’ Also, it is not acting in conjunc- 
tion with a pulse which has reached the A and B synapses from the 
input via the D. neurons (being class (1)); therefore, as seen in 
Section 4, the spectrum of pulses sent to the collector neuron 
simply constitutes ‘‘noise’’ and is suppressed. 

Subclass (1b): This set P,, is significant to the store. Unless 
the ‘‘memory’’ among the thresholds (see Sec. 5, 6, and 7) has de- 
cayed to an insignificant level, this set P,, will act as a decoding 
set for those parts of the message that followed it on the previous 
occasion. 

Class (2): This also is divided into the same subclasses (a) and 
(b). 

Subclass (2a): As shown in Sections 5, 6, and 7, this set P,, in 
conjunction with the input pulse changes a proportion of the thresh- 
olds, leaving a ‘‘memory’’ of it so that on a future reinsertion of 
the set by itself, it acts as a decoding set for the rest of the mes- 
age. 

Sub-class (2b): If the set P,, is still part of the old message 
which is being reinserted, it simply reinforces the memory of the 
old message. If set P,, is now part of a new message then it will 
possibly become linked with the new message rather than the old, 
i.e., in future act as a decoding set for the remainder of the new 
message. Combinations such as this may indeed lead to the ‘‘con- 
fusing,’’ described in Section 9, of one piece of learning with an- 
other similar piece. 

Therefore, for the analysis of ‘‘noise’’ (Sec. 4) we must choose 
a set P in class (1a), while for the analysis of ‘‘memory’’ (Sec. 5, 
6, 7) we must choose a set P in class (Ja) (Sec. 5, case (a)) and 
then a set P in class (2a) (Sec. 5, case (b)). 


8. Consideration of deviations from average. Although this in- 
vestigation has been followed through for average numbers, the 
actual number of pulses in any message length of n pulse lengths 
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will vary from 0 to n since the messages considered consist of 
pulses and pulse gaps of equal probability of occurrence, Thus 
the actual number of pulses in a length of n pulse lengths is sub- 
ject to a binomial distribution deviation. This distribution must 
be considered when fixing 4, the number of pulses that will fire 
neuron G and give an output pulse. So far, A has been chosen so 
that the average noise level does not ‘‘pass’’ the gate G. This is 
equivalent to allowing only information—the recall of a memory— 
and the statistical fluctuations of noise above the average noise 
level to fire cell G. A statistical fluctuation below the noise level 
does not fire cell G and therefore need not concern us. 

Any above-average fluctuation gives a spurious recognition ex- 
cept when it coincides with a desired recognition. Thus A must be 
set so high that the message recalled has only a certain percentage 
of such errors in it, in order that (i) when fed back it is able to re- 
call the next part and (ii) any ‘‘observer’’ can make sense of the 
recall from context in. spite of the errors. 

It is obvious then that provided Q, the number of feedback paths, 
is large enough, the percentage of errors can be made as small as 
desired. But @Q is not necessarily n, the number of pulse lengths 
of message feedback, for the model can be modified by arborizing 
the feedback branches F,, Fy,++++, F _, (Fig. 3) so that they feed 
many more unit stores. In fact, these additional unit stores should 
arise from the generalization of the feedback condition. In Section 
3 the model for convenience was designed so that it contained 
(n - 1) unit stores, the rth unit store having pulse lengths r, r +1 
fed back to it to meet the (n + 1)th pulse length arriving from the 
Dth neuron, Let us suppose now that the model is modified by 
adding additional unit stores so that finally there is in the model 
one unit store for each case r, r + j, where the (r, r + /)th unit store 
has the rth and (r +/)th pulse lengths fed back to it to meet the 
(n + 1)th pulse length arriving from the appropriate D neuron, and 
where 


ALP Oy ives eps and ee at fe eee iy Merida 


Then the total number of unit stores necessary in the model for 
a complete array with a maximum feedback length of 7n pulse 
lengths is N, where 

n(n-1) n?® 


= 


N sn 
5 5 


for n large. 
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The analysis of Sections 4, 5, 6, 7 holds with N replacing (n - 1) 
in the expressions. 
Thus f is given by 


3 N N 
eS os) eee 
2 22, +1 2¢, +1 


Now a measure of the life U of an unreinforced memory in the 
store was found to be given by 


U ~ (22, + 1)?. 


The function U is also a measure of the minimum storage ca- 
pacity of the model since it is the maximum length of message 
that can be stored, no (n +1) pulse lengths of which are repeated. 
It is of interest to obtain some numerical values for the parameters 
involved. For example if 2.10° unit stores are present in the model 


and cees i? the average number of unit stores in each threshold 
@, + 
) 
state,is 107, then U is 4.10° and n = 2.10°. Thus a message could 
be recalled if only 5.107 ° of it were presented as an unlocking key. 


Table 1 gives other examples. 


TABLE 1 


length of size of average number of minimum 

feedback store unit stores in each capacity 

n N threshold state of store 
N/2e, +1 


U 
0 


2.102 2.104 10! 4.106 
10? 4.104 
2.103 2.106 10! 4.1010 
102 4.108 
103 4.108 


A further consequence of deviations from average is that the 
number of limiters changing during a unit of the quantized time 
will vary so that the initial strength of the memory of the pulse 
will itself have a distribution depending upon the relative numbers 
of pulses and pulse gaps in the preceding n pulse lengths of the 
message. Again, when n is large the probability of a deviation far 
from average becomes small. The analysis leading to the expres- 
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sion U ~ (22, + 1)?, representing the life of a memory, is essen- 
tially mmaffected by the findings of the present section. 


9. Properties of the model. Enough is known now about the in- 
formation store to discuss its properties in relation to the basic 
properties of learning associated with living organisms. 

Firstly, it is built round the four postulates listed in the intro- 
duction, viz: 

(a) Physical traces of information fed into the model are left by 
the lowering of thresholds. 

(b) These traces are distributed widely throughout the model. 

(c) The store is randomly connected and need not be subject to 
rigid geometry. 

(d) The information entering the store is distributed uniquely 
throughout the store by the information that immediately preceded 
it and can be obtained at some future epoch by reinserting in the 
store that information that slightly preceded it. 

The first three as we have seen were suggested by the findings 
of neurophysiology. The fourth was suggested by the fact that in- 
formation to be filed requires a filing system, and an inexhaustible 
filing system is available in the information that slightly precedes 
the information to be filed. In addition, a notable feature of learn- 
ing and memory is its chain structure, the recall tending to take 
place in the same order as the presentation of the information. 

In particular, postulates (c) and (d) in their action highlight a 
property that may be characteristic of the storing and recollection 
of information by living organisms. In the present model the stream 
of information alters the internal states of the thresholds of syn- 
apses A and B. This is done in a random fashion since at any 
instant if an observer took the unit stores and inspected their A 
and B synapses, there would be a random distribution of threshold 
states among them. In other words, the entropy of the set of states 
is kept at a maximum. 

Nevertheless it has been seen that a set of pulse lengths which 
has, in feedback, acted as an encoding set (i.e., altered the in- 
ternal states of the thresholds) will be capable, on being reinserted 
in the future as a feedback, of selecting from the ‘‘microstates’’ of 
the maximum entropy distribution a subset that is not randomly 
distributed with respect to the set, and is meaningful. 

The following additional properties belong to it. 
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It is evident that on being presented with any section of a mes- 
sage stored previously, it will recall the rest of the message. This 
is equivalent to reciting poetry, given the title of the poem, or be- 
ing given a couple lines of a well-known poem, recognizing the 
quotation, and completing the poem. 

The memory of a piece of information, unless reinforced at some 
subsequent time, will suffer a gradual near-exponential decay until 
it can be no longer recalled. This decay of memory is one of the 
most fundamental features of learning. Studies show that it pro- 
ceeds most rapidly immediately following learning, then after the 
first few hours the rate of forgetting becomes slower (Miller, 1951). 
This feature is also shown by the present model. 

The store has the property, however, that memories can be re- 
inforced by repetition so that they will not slip below the thresh- 
old of recall. Thus those pieces of information in constant use re- 
main fresh in the memory. Since all practical languages contain 
words and phrases in constant use, the store would build up a 
vocabulary. If a set of phrases were fed into the store incom- 
pletely, in the sense that some of the pulses in the phrases were 
missing, the store would be able to complete the phrases by con- 
text and supply the most-likely following words. This bears some 
resemblance to the ability of human beings to recognize more suc- 
cessfully words accompanied by various noise levels where the 
words are arranged in sentences than where they are rearranged in 
random order and presented accompanied by the same noise levels. 
Work on these lines by G. A. Miller, G. A. Heise and W. Lichten 
(1951) and later by D. J. Bruce (1956) has shown that context in- 
creases the ability of the subject to select the correct word from 
the material presented to him. This ability of context is seen to 
be an automatic function of the present model. 

With the present model, as with human beings (Miller, 1951), 
nonsense messages are more difficult to learn and are forgotten 
more quickly, where a nonsense message is one, the component 
parts of which have been presented once only, while a sense mes- 
sage is one, the component parts of which are in constant use. The 
store behaves in this way since the n pulse length feedback of a 
sense message contains some parts that are in everyday use by 
the store so that the sets of A, B synapses governed by them will 
have already been placed in the one-pulse threshold region. Thus 
part of the work of impressing this message will have already been 
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accomplished. Similarly, these sections are individually being 
often reinforced by everyday use, so that the sense message is 
partially reinforced at intervals and decays more slowly in con- 
sequence. 

The situation, however, is rather more complex than this for in 
a sense the store will tend to forget a sense message by confusing 
it with any other sense message that resembles it, that-is, where 
the second message contains large sections highly similar to sec- 
tions in the first. If the section lengths exceed n pulse lengths, 
the model in recall will tend to follow on such sections with the 
succeeding message most recently presented after that part. This 
is reminiscent of some of the phenomena of retroactive inhibition 
where maximum interference of new learning with old seems to 
occur when the new learning resembles the old (Robinson, 1920; 
Skaggs, 1925). 

One minor feature of the present model is that if the end of a 
message is followed by its beginning without pause, the end be- 
comes the unlocking key for the beginning so that the store will 
cycle that message endlessly through itself in a manner not un- 
familiar to anyone who has had a ‘‘tune running round his head.”’ 
To break such a cycle one consciously thinks perhaps of some 
other tune; to break such a cycle in the model store one would 
have to momentarily raise cell G’s threshold above signal level or 
feed in a new message to change the set of limiter positions. 

If the store was damaged to the extent of having numbers of its 
unit stores destroyed, the effect would be a decrease in efficiency 
in the store’s operations rather than any sudden and complete 
breakdown. Memories in general would not be destroyed. With a 
smaller number of active units, fluctuations from average noise 
level (which would be lower) would be wider and a piece of in- 
formation might have to be reinforced more often before learning 
took place (that is, before a sufficient number of the A, B syn- 
apses that remained were put in the one-pulse threshold region for 
the overcoming of threshold @). Similar results would follow the 
severing of percentages of the feedback fibres instead of the de- 
struction of numbers of the unit stores. If the original number of 
unit stores was large, an appreciable number of unit stores might 
have to be destroyed before a decrease in efficiency was noticed. 
The analogy between damage of this sort and the effect of de- 
stroying or cutting sections of the cortex of the rat is obvious 
(Lashley, 1926, 1929). 


METHOD OF STORING INFORMATION 167 


In Section 7, an expression for the life U of an un-reinforced 
memory in the store was found. The function was also a measure 
of the minimum storage capacity of the model, and the consider- 
ations leading to the table in Section 8 showed that this method of 
storing information gives a large storage capacity. 

The fact that this simple model possesses a number of the basic 
properties associated with learning by living organisms therefore 
seems promising, and it is hoped to continue the study of some of 
the parameters of the model in future papers. 

Finally, something may be said on the unit store itself, in par- 
ticular the limiter. In the simple model used in this paper, two 
Synapses only, connected by a limiter, are postulated. In nature, 
many synapses in general may affect a cell and it would be possi- 
ble to construct more complicated models of the present form where 
the limiter for a many-synapse unit would be a statistical distribu- 
tion of threshold states. A limiter may be no more in fact than the 
principle that, in a cell, where there is only a certain quantity of 
the ‘‘stuff-that-lowers-thresholds,’’ any lowering of one results in 
a raising of the other by impoverishment. The present author is 
not competent to compare further his mathematical constructs with 
living matter. 


In conclusion, he wishes to express his sincere thanks to his 
colleague, Dr. M. W. Ovenden, for many valuable discussions on 
the subject matter of this paper, and to Professor N. Rashevsky of 
the University of Chicago and Dr. E. Trucco of the Argonne Na- 
tional Laboratory for much help and advice in the presentation of 
this paper. 
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We give the solution of a simple linear system of differential equa- 
tions which plays an important role in a recent paper by A. E. Roy (1960). 
The equations are: 


dy, 
rar eet 


dy 


Fe My ~ yt Vy41 (for y = 2, 3,-5.5 m—1) 


dy, 
dt | ¥m-1~ ¥m- 


Their general solution is found to be 
m—1 
kn —B,T 
vy(T)= >) { cos le-» Blas k } 
k=wO 


where 6B, = 4 sin? (k7/2m), and the quantities B, are arbitrary constants 
which may also be written as 


1 m 
By=— ee if 
T=1 


m 
a k 

B,== {008 for-» |} fork =1, 2, heer e Ls 
Tul 


with Y, = y,(0). 
The solution is then specialized to fit the initial conditions chosen by 


A. E. Roy. 


*This work was performed under the auspices of the U. S. Atomic 
Energy Commission. 
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In a recent paper by A. E. Roy (1960) the following problem 
arose in connection with a system of linear differential equations: 
to find the latent roots of the m x m matrix 


-1 1 
lont=s Care. 
1 -2 1 
pee Sh ‘ye | 
(1) 
p oe, a Sa | 
18.1 
Lesage. 
1 -1 
We assume that 
m = 3. (2) 


In Roy’s paper (loc. cit.) the number m was always odd, namely 
m=2e, +1, (3) 


& being a positive integer. We shall, however, not impose that 
restriction, so that m may be an even or odd integer. 

In (1) the elements of the two diagonals immediately above and 
below the main diagonal are all equal to 1; those in the main 
diagonal itself are equal to -2, except for the first and last of 
them which have the value —1; the remaining elements of the 
matrix (1) are zero. 

The characteristic roots, A,, Of (1) are all real, since the matrix 
itself is real and symmetric. Moreover, we must have 


This follows from a well known theorem on determinants which 
states that for any matrix lo, Hl the characteristic roots, con- 
sidered as points in the complex plane, lie within or on the bound- 
ary of the region formed by circles with centers at a;; and corre- 
m 

sponding radii equal to >: a, | (Taussky, 1949). 

j=t 

(i#i) 

The characteristic equation 


A(A) =0 
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is obtained by adding —A to each diagonal element and setting the 
resulting determinant equal to zero. Since the sum of all elements 
in each row of the matrix (1) is zero, it follows that 


A=0 (5) 
will be one of the latent roots. 
Thus, factoring out the term —A, and introducing the new variable 


p= A -2, (6) 


we are left with the following reduced characteristic equation: 


D,(y) = he = 0. (7) 


pH pe Re * 


Let E,, (uy) be the determinant obtained from D,, (y) on replacing 
the last diagonal term, namely p +1, by p Expanding E,, (p) by 
elements of the last row we find that 


Ey (pt) = BE m—1(H) — E mag (H) + (-1)"*". (8) 
A similar expansion for D,, (uy) gives 
Do(p) = (n+ 1) Fn.) - Enea) HC)" (9) 
Equations (8) and (9) imply that 
D,, (pw) = Eg (py) + Em 1 (4) (10) 
and that D,, (uw) must be a solution of the difference equation in m: 
Diy (H) — HD mos (H) + Pm—o(H) = 0- (11) 
Direct calculations for m equal to 3 and 4 show that 
Ds (yu) = Mote 


and - 
D(p) =p - 2b 
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Thus equation (11) holds for m 2 3 if we define 
D,(u) =1 (12) 


Do (py) = bh. (138) 
The general solution of (11), obtained by standard methods, is 


m m 
AC? + Bae; 


and 


where A and B are arbitrary constants while a, and %, are the 
roots of the quadratic equation 


a?-pa+1=0, 
1.85: 


(14) 


with 
amps 4. (15) 


The solution satisfying the initial conditions (12) and (13) is then 
found to be 


D(H) = = c + 8)" —(n- 2" 


(16) 
(m = 7 3S, ees hs 
We notice that 
D, (—H) = (-1)"*? D,, (p) (17) 
and 
D (£2) = m-(£1)"71. (18) 


The right-hand side of equation (16), which is a polynomial in yp, 

n be transformed as follows. 

Expanding (y+ s)™ and (1 — s)™ by means of the binomial theo- 
rem, and subtracting term by term, we obtain 


af x m 
D(z) ae EL > bai ‘ :) p™-1-2p 32h, (19) 
P=0 


where «x is the largest integer smaller than Ma If now we use again 
2 


the binomial theorem to expand 8°? = (y?~ 4)? in powers of u 
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we find, after rearranging the resulting double sum, that 


1 K 
D(z) = —s (-1)” Q2v C,, edt pat (20) 


v=0 


with 


Se m pel AE m 
3, () (pe) Daal o hs iereee) py 


But the right-hand side of equation (21) is equal to 


gm—1-2y oP mt as 
V ] 


as follows from equation (31) on page 253 of E. Netto’s Lehrbuch 
der Combinatorik (1901), after some appropriate changes in notation. 
We thus obtain the final result 


D - vy {[m-1i-v 
nw) =) 2) ( ymoinay 
V=0 ‘4 


m—1 


if m is odd (22) 
with = assy eee 
3 —1if mis even. 


For the purpose of finding the latent roots, however, it is more 
convenient to start from equation (16). As we have seen, we need 
only consider real values of y such that 


9 2 p< 42 (23) 


(cf. (4), (6), and (18))- 

If (23) holds, the value of 's given by (15) is purely imaginary. 
We now introduce a quantity 9, defined by 
(24) 


cos 9 = -, 


wolr 


0<pK<z. (25) 


Then 
pts=2 ete. (26) 
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so that equation (16) becomes 


sin (mq) 


D,,(u) = —- (27) 
Sin © 
D,, (u) vanishes only when 
sin (mp) = 0. (28) 
Thus, by (28), (24), and (25), the m-1 roots of equation (7) are: 
7 
py = 2 coss=— 
ze (29) 
with fen Set sae RE egos me 
If we put 


as was done by Roy in Joc. cit., we obtain for the latent roots of 
the original matrix (1): 


k r 
B, =2 [1 — cos =) = 4 sin? ees 
m 2m (31) 
with f= 055159. 05g eae 
The smallest of the positive B’s, namely B,, is given by 
Tv 
B,=2 [1 ~ 0s =|. (32) 
If mis sufficiently large: 
7 
aT at ea Ry eps 
m 2m?’ 


and equation (32) becomes 
B, a aS * (33) 


We may now proceed to solve the system of differential equa- 
tions with which Roy was concerned. Using a notation slightly 
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different from Roy’s it can be written in the form 


dy, 

dt 77M +% 

dy, 

ae = yoy SY + Vyy (for v =Q; 8,....55:m—1) (34) 
dy 


As is well known, we have to find a non-trivial solution 
Gis Meineke oy 
of the equations 
(u, +1) a4, +4, =0 
Gy spt HE SyE t+ Fy, = 9 (vy = 2; 3,..<, m—1) (35) 
Ps ehaet 2) Sek 0 


for each value of k(k=0,1,........ m—1); here the p;’s are the 
quantities given by equation (29) and, in addition, py, =-2. The 
solution of (35) will be uniquely determined once a,, is given. 
For then we can calculate a,, from the first equation, a,, from 
the second equation, and thus all the a,,’s follow successively. 
Choosing 


kn 
ask = COS om’ 


one can check without difficulty that the quantities a,,, as ob- 
tained from the system (35), will be 


kn 
ay, = COs (2 v -1) am 

for ped Oot Gm (36) 
k=O) 4h isi. m1; 


Hence the general solution of the differential equations (34) is 
given by 


m-1 
¥,(T) = ne Bows, e7 Pkt } (v=1, 2, cong m), (37) 
k20 
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where the @,,’s must be inserted from equation (36), the B,’s from 
equation (31), and the quantities B, are arbitrary constants, to be 
determined by the initial conditions. It can easily be verified 
directly that the expressions (87) satisfy the system (34). 

Finally, we may wish to express the constants B, in terms of 
the initial conditions. Let 


Vy (Oe 1) (ele «le (38) 
We then have to solve the linear system 


m—1 


a,, B, = ¥. (ede? pa. ty. m1) (39) 


Vv 
k=0 


for the unknown quantities B,; we recall that the coefficients Qi, 
are given by equation (36). 

Since the roots of the characteristic equations were all real and 
distinct, the m sets of functions, 
—pat —h7% - tT 
e Tn et at : yeh) « Seguin a Bm—1 


(ib, a. one, ite 


a0 


are linearly independent. This means, in particular, that the de- 
terminant la. must be # 0 (see, e.g., Kamke, 1956, chap. Visite 
verify this conclusion directly, let us call D, (c) the oxo de- 
terminant 


D (0) = det |cos {im +-o)-11 Ar} 
2m (40) 
with v=ul, 2, a5 o and £=0,1, .:.5 6 —4, 


o being an integer such that 1<o<m. We have D,,(1) =1, and 
the determinant of interest is D,, (m). 

For each column of D, (co), and for p = 1,2,....0-2, add to the 
vth element the (v +2)nd element, as well as the (v + 1)st element 


multiplied by 2 cos C +a-1) |; to the (o —1)st element add 


the oth, multiplied by 1 +2 cos G +oa-1) |. These trans- 
m 
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formations do not change the value of D, (c), but they lead to the 
recurrence relation 


D(a) = (2173 .920-3 .sin |(o _ 1) =| ° 
m 


20=—3 ‘ 
Ay eee i (41) 
(a ee fe, Dm(o-1). 


The factor multiplying D,,(o-1) in equation (41) is not zero if 
20m, and so 


D,,(m) #0, (42) 


as expected. 

The equations (39) therefore have a unique solution; to find it 
we must obtain the inverse of the matrix lle,, |, i.e. we have to 
determine quantities 6, such that 


D>, Gun bur = Oye 
k=0 (43) 
for o_o Gah ie ere 


where 
St Ul oy =f, and 6,,=0: if wr. 


We now assert that the elements of >, | have the following 
values: 


and (44) 


Kee 1 Bo aso pthread, 


2 kr 
bia eT al? 
a Se ler va he een ean 


To verify that the righi-hand sides of (44) satisfy the relations 
(43) we may employ the well-known equation 
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cos 9+ cos(@ + 6) +cos(6+25)+...+cos[6+(n-1)8] 


= Sean eA E + eS AD ek (45) 
5 ye 


(For proof see, e.g., Siddons and Hughes, 1928, pp. 263-4). 
Taking for a,, and b,, the values given by equations (36) and 
(44) we obtain 


A Migore sk tk 
Da OO, = ee yi {eos ee Pee a (46) 
va PI m m 


S=vir-—1, 
(47) 
t=v—-7 
The integers s and ¢ are such that s is odd when ¢ is even, and 


vice versa. Using this fact, as well as the relation (45) with 
appropriate values of 6, 5, and n, equation (46) finally reduces to 


with 


—— if wee 
m-1 m 
pk xayer (48) 
km m—1 
ifv=r 


On adding a,b, = J to both sides, equation (48) becomes identi- 
m 


cal with (43), which we wanted to prove. 
The solution of our system (39) can therefore be written 


By =). o.7, (49) 


with the quantities b,, given by equation (44). We notice that 
Peymaimr = (~1)* by, (50) 


Hence the constants B,, with even k depend only on the sums Y + 
AE whereas the B,’s with odd index k are linear functions of 
the differences Y-=Y 


mt l=r* 
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Equation (37), together with (31), (36), (49), (38), and (44), gives 
the complete solution of the system (34), expressed in terms of the 
initial conditions. 

Let us now consider the special initial values chosen by Roy in 
loc. cit., namely 


¥, =f) 

2 

Y= 3 So (iim Bad stm ~ 1) (51) 
1 

Tn = 3 Sos 


where S, is a constant. 
It is found that all the quantities B, with even & vanish, except 
for B, which is equal to (2/3)S,. On the other hand, for odd & 
4$ k 
B= =— cos alt 
3m 2m 
We assume at this point that m is an odd integer (cf. equation ; 
(3)); after suitable transformations the solution for the particular 
initial conditions (51) becomes 


m= 1 
2 : ey pee 
¥y(T) = Bo {2 . x Cyje 3 fp (52) 
j=1 
where 
of Ae +3 
Cy; = ee [eos Cintue Peon Bieta 
“A (53) 
4 [ei Mere) (Q}—1)9 
Be EIN ei bee Sg aii Bol PS Pens heel has 
3m 2m oo 
and 
ae 
2 —_—_- CC" 
B,=4 sin I ee : (54) 
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A model of the processes occuring in the exchange of a drug between 
capillary plesma, extracellular space and intracellular space is devel- 
oped. This leads to an interesting set of differential difference equa- 
tions, one of which is an integrodifferential equation, another a partial 
differential equation. Under certain conditions, these may be simplified 
to a set of ordinary differential equations. The application of Laplace 
transform techniques to the solution of these equations is discussed, 


1. Introduction. Mathematics plays a fundamental role in the phys- 
ical sciences because it furnishes the possibility of setting up 
realistic models of physical phenomena which can then be treated 
by uniform techniques. 

This has not generally been true of the fields of biology and 
medicine. Often the realistic models have defied the state of the 
art of mathematics, or the oversimplified models that have been 
treated have furnished little insight into actual biological phenomena. 

The present paper is the first of a series of papers in which we 
wish to study qualitative and quantitative aspects of chemotherapy. 
Specifically, what we wish to study is the distribution of a com- 
pound in the organs of the body after its injection into the blood 
stream. 

These problems are similar in many respects to a class of prob- 
lems arising in mathematical economics in connection with mul- 
tistage production processes and generally in the domain of ‘‘in- 
put-output’’ analysis. The biological problems are, however, in 
many ways more complex. 
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It turns out that the study of either dynamic or static (steady 
state) behavior of the organs of the body, arranged in series-par- 
allel networks, bears a certain resemblance to the operation of 
complexes of interdependent industries, and also to the cascade 
processes encountered in isotope separation or in the operation of 
oil refineries. These matters will be discussed subsequently. Here, 
in this initial paper, we shall study a relatively simple system con- 
sisting of the heart and one organ. 

At the moment, we are concerned with the purely descriptive 
phase of the study of these life processes. As we shall see, one 
of the advantages derived from the formulation of mathematical 
models on even this moderate plateau of realism is the light it throws 
on the need for further experimental work which will yield the data 
required for actual numerical results. 

At the moment, we wish to consider only the analytic aspects. 
Subsequently, we wish to formulate multi-organ models and apply 
analogue and digital computing techniques to the solution of var- 
ious linear and nonlinear sets of equations which arise from them. 

Following this phase of our program, we wish to turn to our pri- 
mary objective, the study of control aspects of these processes. 
Our ultimate aim is to determine chemical injections with optimal 
specific properties, 

From the mathematical point of view, the study of these prob 
lems is quite interesting due to the presence of novel features not 
usually found in the study of mathematical physics or mathematical 
economics. To begin with, there are a number of difficulties pres- 
ent at the very outset when we attempt to set up a mathematical 
model. Next, we find that the built-in periodicity and time lag fur- 
nished by the circulation of the blood leads to differential difference 
equations involving partial derivatives, Furthermore, the time lag 
appears in the boundary conditions! Although these can be treated 
by Laplace transform techniques under certain simplifying assump- 
tions which yield linear equations, more detailed models yield non- 
linear equations. The numerical solution of large-scale systems 
of differential difference equations is not a trivial matter, In par- 
ticular, the computer memory requirements may become excessive. 
The equations resulting from the Laplace transformation are also 
of interesting type. They are integrodifferential equations with 
boundary conditions of two-point type. 
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The question of determining what type of steady state results 
hinges upon the investigation of the characteristic roots of quite 
complicated transcendental functions. Although the method of L. 
Pontrjagin (1955) can be used in some cases, in others, the matter 
is still unresolved. As we shall see, physical reasoning can some- 
times be employed. 

The reader versed in biological and physical matters will see 
that we have omitted a number of effects in our treatment, Our ex- 
cuse is a simple one—the problem is a complicated one even at 
this level. 

In Section 2 some physiological considerations are provided by 
way of background. This leads to a model of a one-organ being in 
which we seek to determine the concentration of a drug as a func- 
tion of time and position, where the concentration satisfies a cer- 
tain partial differential equation along with appropriate initial and 
boundary conditions. The role of the Laplace transform in effecting 
a solution is indicated. It is then shown in Sections 9 to 11 how 
many essential features of this original model can be preserved, 
though the problem is much simplified, through the introduction of 
a still further simplified model. Finally in Section 12, we show 
how to formulate equations for the situation in which a drug is in- 
jected into the blood stream, enters the extracellular space, and 
finally enters the intracellular space where it takes part in a re 
versible chemical reaction. 

2. Physiological considerations. Except for two special cases, 
we may consider the various organs of the body as being in contact 
with the closed-circuit circulation as in Figure 1, with the various 
organs or regions of the body in parallel. 

As an introduction to the quantitative study of this system, we 
consider the drastically simplified case of a single region with a 
closed-circuit blood supply. We defer consideration of mixing in 
the circulation until we take up the many-organ model (Jacquez, 
Bellman, and Kalaba, 1959). We feel that this case will illustrate 
the essential assumptions required to construct the model. 

The region involved consists of many simply connected spaces, 
living cells, of various roughly similar shapes and sizes; each is 
separated by its boundary, the cell membrane, from a watery fluid, 
the extracellular fluid, which fills the space outside the cells. The 
blood supply to this region enters via a number of tubes, arteries, 
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HEART 


ORGAN 3 


FIGURE 1. Diagram of parallel arrangement of body organs in relation 
to blood supply. 


which branch into the small arterioles. Each of the many arterioles 
gives rise to many, small, thin-walled tubules, the capillaries, 
which pass in many directions through the extracellular space and 
finally join into venules which combine to form the veins which 
carry the blood from the region. Exchanges between the blood and 
and the cells of the region occur by diffusion between the blood and 
extracellular space and between extracellular space and inéracel- 
lular space. 

Although capillaries vary somewhat in dimensions, we may take 
the average capillary to be .0004 cm in radius and .04 em long. All 
the capillaries to a region are not always open with the number 
through which blood is flowing at any one time depending on the 
state of the tissue. Thus, for a resting muscle only 5 capillaries 
per sq. mm may be open, while for stimulated muscle as many as 
190 capillaries per sq. mm may be functioning. We consider only 
situations in which the number of capillaries open is roughly con- 
stant. The blood flowing through the capillaries consists of a sus- 
pension of cells in fluid, the plasma. The exchanges we will be 
concerned with will be between the plasma and extracellular space, 
So we will speak in terms of the plasma flow to a region. 

We may now construct our model. The fully rigorous approach 
would be to write the differential equations of material balance 
taking into account the flood flow rate for each capillary, filtration 
and diffusion across capillary walls and diffusion in the extracel- 
lular space, and introducing the geometrical relations between al] 
capillaries in the region and the distribution of extracellular space 
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in relation to cells and capillaries. The explicit statement of such 
a program is sufficient to indicate its sterility. A macroscopic 
viewpoint may be more useful. If one takes a cross section of a 
region, the capillaries in the cross section will be roughly uniformly 
distributed with their flows in various directions. From a macro- 
scopic viewpoint, such a distribution tends to smooth out any dif- 
fusion gradients. The normal rotation of functions among various 
capillaries so that capillaries without any flow open up while others 
shut down can only add to this leveling process. Filtration across 
capillary walls may be neglected since it is a minor process com- 
pared to diffusion, for molecules smaller in size than inulin (Chinard, 
Vosburgh, and Enns, 1955; Pappenheimer, Renkin, and Borrero, 
1951). The flow rates for various capillaries are not necessarily 
the same. We do not introduce this in the present treatment but 
assume one average flow rate. These considerations suggest the 
most important assumption of the model, that is, at any time the 
concentration of a material exchanging between the plasma and the 
region is uniform in the extracellular space. The approximate va- 
lidity of this assumption is implied by the experimental evidence 
that the tissue-blood exchange rates for inert gases (Jones, 1950) 
and small ions (Conn and Robertson, 1955; Dobson and Warner, 
1957; Freis, Higgins, and Morowitz, 1953) are limited by blood flow 
to a region and not by molecular diffusion. Thus we may concep- 
tually lump all of the capillaries into one the length of which is 7, 
the average length of a capillary. As with the flow rates, we dis- 
regard the effect of any possible distribution in the lengths of cap- 
illaries. The effect of this will be considered in a future paper. 
As shown in Figure 2, this capillary contains a volume &, of plasma, 
has a surface area 4, of contact with the well-stirred extracellular 
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FIGURE 2. Model of a closed circulation with a single organ. 
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space of volume #,, and has a volume flow rate of c. The extracel- 
lular space is in contact with the intracellular space of volume fF; 
and surface area A,. . 

G. W. Schmidt (1952, 1953) has considered a model similar in 
some respects to ours for the exchange of a substance between 
blood and extracellular space, but has not explicitly introduced 
the problem of recirculation or of the possible reaction of the sub- 
stance with a component of the intracellular space. R. E. Smith 
and M. F. Morales (1944) and Morales and Smith (1948) have treated 
the problem of the exchange of an inert gas between capillaries and 
tissues. W. C. Sangren and C. W. Sheppard (1953) have presented 
calculations on a model which assumed rapid mixing in the direction 
perpendicular to the axis of the capillary but no mixing longitudi- 
nally. Considering the anatomical distribution of capillaries and 
extracellular space this does not seem to be a physically reasonable 
assumption. 

We consider first the simplest case, that of diffusion of a com- 
pound into the extracellular space. The compound is assumed not 
to enter the intracellular space. We shall subsequently consider 
both diffusion processes, as well as an active transport process. 


3. Notation. To reduce some of the foregoing ideas to mathemati- 
cal form, let us introduce the following definitions: 
u(@, t)—concentration of drug in mass per unit volume at 2 in the 
capillary at the time ¢ 
v(t)—spatially uniform concentration of drug at time ¢ in the 
extracellular space 
c—volume rate of plasma flow in the capillary bed 
k—permeability constant for capillary walls in units of J¢7! 
J—length of capillary 
f,—plasma volume in capillary 
f&,—volume of extracellular space 
A,—total surface area for diffusion between plasma and extra- 
cellular fluid 
T—the time for a particle to travel around the circulation from 
1 to 0. If R is the plasma volume other than that in the 
capillary, ct=R. We neglect mixing in the heart and 
large vessels. This will be considered in a subsequent 
paper. 
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4. Derivation of basic equations. Let us now derive equations 
connecting u(z, ¢) and v(t). From conservation principles, we ob- 
tain the relation 


R R 
Foe t+h)dz = eat (a t)dx + (u(x, t) —u(e + dz, t)) he 


(1) 


(t) — u(a, t)) hdz. 


The second term on the right-hand side arises through the motion 
of the blood in the circulation and the last through diffusion into 
the extracellular space. Letting A and dz tend to zero, we obtain 
the equation 


R 


Pp 


U 


where u, and u, denote partial derivatives with respect to ¢ and @. 
exominine the rate of change of the compound in the extracellu- 
lar space, we obtain the relation 


kA,A f! 
Riv(t +h) =R,v(t) + e | [u(z, t) - v(t)] de. (3) 
) 


The integral represents the net gain in the drug in the extracellular 
space during the time 4. This leads to the equation 


kA I 
Ro (t) = : 2 | u(a, t)de —kA,v(t), (4) 
0 


valid fort >0,0<2#<l. 
To approximate an intravascular injection we assume that a 
rectangular concentration wave enters the capillary at ¢ = 0. 
Thus, we have for initial conditions, 
u(@, 0) = v(0) = 0; (5) 
while for boundary conditions, we have 
to, OxtSt, 
u(0, t)= 40, t,<#S57 (6) 
u(l,¢-T), LT, 


where o/c <t,<T. 
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The restrictions on ¢, are introduced in order to approximate an 
injection of duration less than one circulation time. The assump- 
tions involved in the derivation of the equation are probably less 
justified if we attempt to describe events of very short duration, 
particularly at the beginning of the process. This would be true of 
an injection of such short duration that the volume of plasma oc- 
cupied by drug represented a small fraction of the volume of the 
capillary bed. For this reason the volume of plasma occupied by 
drug (c¢,) is made larger than the volume of the capillary bed, (R,)- 


5. Discussion. The problem of equations (2), (4), (5), and (6) is 
of an interesting and rather unique type. In particular is this true 
of the boundary condition, involving as it does a time difference. 

To resolve this system of equations, we shall proceed formally 
under the assumption of existence and uniqueness of solution. Our 
principal tool will be that general factotum of analysis, the Laplace 
transform. 


6. Solution by Laplace trans form—I. To simplify the notation, let 


a 
i t (7) 
= ct. (8) 
Substituting in (2) and (4), we obtain the equations 

Roug=-u+K(v-w), (9) 

1 
R (8) = ial udy - Ko, (10) 

0 


for @>0,0<y<1. Here K is the dimensionless constant kA, 7OX 
The sitial conditions then become 


u(y, 0) =0 1 
v(0) =0 a 
while the boundary conditions become 
u,9<0<86,, 
u(0, 2) = <0, 0,<O<T, 6, =ct,, fT =e7% (12) 


u(1, 6-1), 2@>T. 
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Introducing the functions 


L [u(y, 6)] = iG u(y, de °°d6 = Uy, 8), 
E (13) 


L [v(@)] = | v (d)e~$°dd = V(s), 


0 
we see that (9) and (10) yield the equations 
RosU=-U +K(V - 0), (14) 


1 
RisV = Kf Uy, s)dy — KV, (15) 
) 


The boundary condition in (12) transforms into 


U(0, s) = = P= at *1).4-6 2720 (4, 8). (16) 


7. Solution by Laplace transform—lI]. Eliminating V, we see that 
U, as a function of y, satisfies the integrodifferential equation 


k? 1 
U sR_+K)U — ——_— U(y, s)dy = 0. i 
ph, +U- Fae | Ua sdy=0, 7) 
with the boundary condition 
Uu = i 
U(0,s)=-2[1-e “*] +e U (1, 8). (18) 
8 


Let us suppress the s-dependence momentarily and consider the 
problem of solving a linear functional equation of the form 


1 
tb(y) + aw(y) = [ w (y)dy, (19) 
0 


w(0)=c + dw(1). (20) 


A rigorous discussion would take us too far astray. Let us assume 
the existence of a unique solution and see how to obtain it, Set 


iD w(y)dy =m, (21) 


10) 
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an unknown constant, Then the solution of 
w(y) + aw(y) = bm (22) 


is given by 
_ bm 
wy) =fe + . (23) 


where f is a constant of integration. To obtain a relation between 
f and m, we return to (21). The result is 


ey Ee ea (24) 
a a 
or 
f( 6°) 
a hae By CG?) 


To determine the remaining constant, f, we turn to the boundary 
condition of (20). Having obtained the value of f in this way, the 
function w(y) is completely determined. 

Upon carrying out this program we find, after considerable com- 


putation, that 
u -s9 
“2 [l-—e 5| 
8 


UY, 8) = —=——$$$_$—_____ 
b Sas ea 
c ¥ Fins z a (26) 


x jen a ame] : 
a(a—b) 


2 


where 


ei 5) a2 BF +Ks(R,+R ) 


a@=sR +K, b= ; 
4 sh +K sR, +K 


(27) 


8. Solution by Laplace transform—lII]. The desired solution 


u(y, #) is obtained by taking the inverse Laplace transform of 
Uy, 8), 


1 
u(y, 0) maa | UG sje ds, (28) 
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where ¢ is a contour lying to the right of all the poles of U(y, 8) 
considered as a function of s. 

One pole will always be at s = 0, corresponding to the steady- 
State solution, It is easy to see on physical grounds, but difficult 
to establish by analytic techniques, that the other poles will have 
negative real parts. The inversion of equation (26) is a difficult 
program analytically. However, we may check the steady-state so- 
lution by evaluating the residue at the pole s = 0, which yields, in 
view of equations (26) and (27), 

Peeeety eyo eg (29) 
s>0 ave R, + et 
The same result can be obtained on physical grounds from a con- 
servation condition as follows: 
ct 
lim u(z,¢) =u, pane eee, (30) 
t 00 Ro + Re bet 

9. A simplified model. In an effort to relieve some of the mathe- 
matical difficulties associated with the previous model, let us 
handle the z-dependence in a simplified form. We note that the ap- 
pearance of the integral in equation (3) implies that diffusion into 
the extracellular space is determined by the difference between the 
spatial average concentration in the capillary and the concentration 
in the extracellular space. We assume there is a certain concen- 
tration of the drug at the arterial end of the capillary and another 
at the venous end. The diffusion into the extracellular space is to 
depend on the difference between the average of the venous and 
arterial end concentrations and the concentration in the extracellu- 
lar space. This again makes the diffusion dependent on a value of 
the plasma concentration intermediate between the concentrations 
at the ends of the capillary. Figure 3 illustrates this situation. 


FIGURE 8. A simplified model. 
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wu, (¢)—concentration of the drug at the arterial end 
u.(t)—concentration of the drug at the venous end 
v (¢)—concentration of the drug in the extracellular space 

The equations which describe the development of the process 
are: 


lio, O <7 S7,, 
tt (Dace One eer (31) 
1 a b 1 = = 

u(t¢-T),¢>T, 


where Rp/c <t,<T. 


R,b(t)=h,A, ee) ~y 0) (32) 
R, pee) =0(u,-u,)- A, OBE we a) . (33) 


Making the substitutions @ = ct and K = kA /e, one obtains 
Uy» D<O< 4 
u,(9)= < 0, a= 0 <7, (34) 
u(@-T), O>T, 


where 6, = ct,, T = ct. 
R, jess Sa =(u, —U,)-K ea -»| (36) 


With the initial conditions 
u,(0) = w,(0) =v (0) = 0. 


The approximation to the original set of equations is poor for 
processes of duration less than the capillary transit time; this is 
particularly true for those periods when the Square wave is enter- 


U, +Ue 


ing or leaving the capillary bed, for then is @ poor approxi- 


l 
mation to vt | u(@, t)dt. This is also a poor approximation for 
0 
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K very large. For K >2 this model leads to negative concentra- 
tions in the plasma leaving the capillary bed for some short period 
as the square wave enters the capillary bed. This may be seen 
from the following argument. 
For ¢ = 0, u,(0) =v (0) =0 
Thus equation (33) reduces to 


Rp io) _ KA 
a (0) = cu, tea (37) 


which rearranges to equation (38): 

(0) = Bes (38) 
Thus if K > 2, (0) <0, and since u,(0) = 0, this would give nega- 
tive values for u,(¢) for some short interval at the start of the proc- 


ess. In the following section we assume K < 2. 


10. Solution by Laplace transform. We introduce the Laplace 
transforms 


L (w,(@] =U,(s), (39) 
L [u,(6)] = U,(s), (40) 
L {v(6)] =V (s). (41) 


In view of the equations of Section 9, these transforms satisfy 
the equations 


Us) = 2 [1 - 2 Pat's e~*7U,(s), (42) 
R,sV (s) =K pina TM) av (s a), (43) 


tat (8) + U,(8)] - 


—2 = (s) - U,(2) ~ ee 


A) ae 2). 


If next we solve the system of equations (42), (43), and (44) using 
Cramer’s rule, we find that the determinant D(s) which occurs in 


(44) 
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the denominator is 


iL aa eg} 0) 
K K 
ae 2 —(K +s.) 
Daa3 2 ( (45) 
K R 
| (1 +2 +2) - 
2 2 
or 
. io ade ay (% 3% nee 
=——+ + 8! —+——- g - 
(3) Pig ee AAG 5 
(46) 


K? ere R) é sk, ) 
SS ee } a : 
5 e + ( + 8k. 5 5 


Since w,, u,, and v are bounded, on physical grounds, D(s) can- 
not have any roots with positive real parts. The origin, of course, 
is a zero, aS one sees by inspection. This gives rise to the con- 
stant term in the inverse transforms. 

An investigation of the precise location of the roots of D(s) 
would provide valuable information concerning the asymptotic be- 
havior of the functions wu,» Uy, and v. We conjecture that the root, 
other than the origin, which lies closest to the imaginary axis in 
the s-plane lies on the negative real axis. If this is so, then the 
functions u,, w,, and v have the form a 1 7 @2€ Raves t 0, 
where a, is the root, and a, and a, are the residues at zero and a, 
respectively, Experimentally ee constants could be determined 
by recording the concentrations as functions of time and then se- 
lecting the constants to provide the best fit to the experimental 
curves, 

The location of the least negative root of D(s) can easily be ac- 
complished graphically. Equation (46) can be rewritten as 


3 KR, KR 
-—RR vs ( ~) o 
ore 8 : 2 
ee RK R,K Sat 
ede ee 9 a 
so that it is merely necessary to graph the exponential curve of the 


left-hand side and the rational function of the right-hand side and 
find their intersection points. One lies at s = 0, and the other is 
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negative, as an investigation of the roots of the polynomials in the 
numerator and denominator of the fraction in equation (47) shows. 
It can also be seen that an increase in T decreases the rate at 
which equilibrium is approached, which is physically quite reason- 
able. 

The location of the roots of exponential polynomials such as oc- 
curs in equation (46) has been the object of many investigations 
(Ansoff and Krumhansl, 1948; Collatz, 1947; Hayes, 1950; Lax, 
1948). In particular, we should like to call attention to the results 
of Pontrjagin (1942), which are described by R. Bellman and J. M. 
Danskin (1954) and are available in translation (Pontrjagin, 1955). 


11. Numerical aspects. The system of equations (31) to (35) 
lends itself well to solution using a high-speed digital computer. 
The function u is known on the interval [0,T], and w and w can be 
determined on this interval using the differential equations (32) and 
(33) along with the initial conditions of equation (35). Then equa- 
tion (31) yields w(¢) for T S$ ¢ S$ 2T, and so on. Having the time lag 
present increases the size of the computer memory required, but is 
otherwise innocuous. For more realistic models, this increased 
demand for memory space could become important. 


12. A model for intracellular penetration. As was stated earlier, 
we wish to consider the case in which the drug enters the cells 
themselves and there reacts chemically with substances within the 
cells to form new compounds. In general, these reactions may be 
reversible. The situation which we wish to discuss is shown 
diagramatically in Figure 4. 

The function u(z, ¢) is the concentration of the drug in the capil- 
lary at position 2 at the time ¢. The functions v(¢) and w(t) repre- 
sent the spatially homogeneous concentrations of the drug in the 
extracellular and intracellular spaces, respectively, and 2(¢) is the 
concentration of the compound ED which is formed from the chemi- 


CAPILLARY 


ARTERY 


FIGURE 4. An intracellular penetration schematic. 


196 R. BELLMAN, J. A. JACQUEZ, AND R. KALABA 


cal union of the enzyme F and the drug D. We assume that both 
the enzyme F and the compound EP do not diffuse out of the intra- 
cellular space. We do wish to include the possibility that there is 
an active transport process which tends to concentrate the drug 
intracellularly. The transport term is assumed to be linear in the 
extracellular concentration. This is only a fair approximation even 
for-low concentrations for the amino acid transport system (Heinz, 
1954; Jacquez, 1957). It is used for the time being to demonstrate 
the effect of an active transport term. 
The equations for this model become 


R kA 
FU ee yeiien Ms (48) 
kA mere 
Rv = 7 / ula, t)dx ~ kA,v(t) - mA(v - w)-nAn, (49) 
0 
fiw = mA(v - w) +nAyv + k,R2—-k Rw [E, - 4], (50) 
4=k,w [k, -2]- kz. (51) 


Equation (48) has been met before. In the next equation, the last 
two terms measure the rate of decrease of v as a result of both 
diffusion and transport into the intracellular space. The last two 
terms in equation (50) measure the rate of increase of the concen- 
tration of the drug in the intracellular space as a result of the de- 
composition 


ED *3, p +E, (52) 


and as a result of the reaction 


E+D “5 ED. (53) 


In these equations, z is the concentration of the compound FD. 
Following the idea of Section 9, a simplified version of this prob- 
lem could be considered, 


13. Discussion. The mathematical equations obtained from the 
model of drug distribution in a one-organ entity present obvious 
analytical difficulties. As will be shown in our next paper, these 
are multiplied many-fold when one attempts to link a number of 
such models via a circulatory system to give a model of a many- 
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organ entity. It is to be expected that modern computational tech- 
niques will help clarify the nature of the solution to such systems. 

However, the purpose of such an enterprise is more than to ob- 
tain an explicit mathematical solution of a model. As with all 
theorizing, it represents an attempt to gain understanding of a com- 
plex process. With this added insight, we then return to experi- 
mentation, and test various hypotheses. The results of this testing 
enable us to build more complex models, leading to further experi- 
mentation. We thus have a complex feedback process. 

The very act of setting up a mathematical model points the way 
to experiments, to the need to measure the parameters which appear 
in the equations: in this case, the blood flow and blood volume of 
an organ, the capillary permeability and area for diffusion in an 
organ, and the permeabilities of various cell types. Thus even the 
formulation of a model which is a crude approximation to reality 
may be useful because it provides a well defined basepoint from 
which to strike out in search of new hypotheses and designs of 
experiments. 
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An error appearing in the proof of Theorem 4 of a previous paper of the 
author’s (1959, Bull. Math. Biophysics, 21, 289-97) is pointed out, and a 
new proof of the theorem is supplied. We also obtain a corollary from 
Theorem 3 of loc. cit. which reveals the existence of a hitherto unrecog- 
nized class of codes. 


In the proof of Theorem 4 of our paper ‘‘Some Further Comments 
on the DNA-Protein Coding Problem’’ (1959b, Bull. Math. Biophys- 
ics, 21, 289-97), familiarity with which is presupposed, we use in 
an essential manner the inequality labeled (6) appearing on page 
294; namely, 


Se rr, (wv) y r Up, + €)A(w)] 
1 


Sse ssa sk Ses Eee (1) 


s (r; Be 8;) d,(w) y (r; at 8) l(p, + €) A(w)] 


It has been pointed out to the author by Mr. Richard White (Depart- 
ment of Physics, University of Wisconsin, Madison) that this in- 
equality is incorrect. In fact, if (1) holds, then we also have, by 
the same argument, 

nt n 
y. 8,4, (wv) s s,[(p, + &) Aw)! 
1 1 


<—. ?) 
>, (r+ 8A, (w) >, + I; + e)A(w)] 
1 1 


199 


200 ROBERT ROSEN 
Adding (1) and (2), we arrive at the absurdity 


a (r, + 8,)A,(w) Sy (7, + 8,)[(p, + €)A(w)I 
[2 
d+ 8), (w) DY" (7, + 8) Up, + &) A(w)I 


Hence the remainder of the argument, which depends on the in- 
equality (1) above, is incorrect, and Theorem 4 is not proved. 

Nevertheless, the statement of Theorem 4, as it appears in loc. 
cit., is correct. We shall now provide a valid proof of this state- 
ment, which may also have some independent interest. We wish to 
thank Mr. White and Professor N. Rashevsky for examining this 
proof, although the author assumes full responsibility for any short- 
comings which may remain. 

Let us denote by F the closed unit interval, and by ?, the carte- 
sian product of & with itself n times. That is, R, is the closed 
unit n-cube, #, is the unit square. We can construct a mapping Q: 
I, — A, by writing 


“ r,(w) A,(w) A, (w) 
ete) (tay atop ane) 
and a mapping ¥: Il, — R, by writing 


u(y) = no(v 
YW (r) = (a0 ’ 2 ), 

u(v) wv) 
where the symbols have the same meaning as in Joc. cit. If now 
f: 0, — 2, is a monomorphism, then we can construct the follow- 
ing diagram: 


Here, f is the unique mapping which makes the diagram commuta- 

tive; more explicitly, we have f(9(w)) = UW (f(w)) for each weMN.. 
We assert that Theorem 4 will be proved if we can show that ii 

can be extended in a unique manner to a continuous mapping, con- 


tinuity being relative to the usual metric topology on RF, and k,. 
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Before we proceed to justify the assertion just enunciated, let 
us clarify its meaning. We first observe that the range of ~ (which 
we shall denote by r(~)) is dense in some subset of F 3 likewise, 
r(w) is dense in a subset of R,. That is, given any Boi p=(p,, 
Pov+++, D,)Er(P) and any €> 0, there exists a point gér(~) such 
that g lies in an €-neighborhood of p. 

To prove this, we note that, according to Theorem 4 (Rosen, 
1959a), pe€r(P) means precisely that 9~'(p) is a submonoid of ais 
in fact, this is the ergodic submonoid that we denoted there by 
Evp;)» Hence, ? —'(p) contains words of arbitrarily | large length. 
If now €>0O is given, we can find a word uw in ~~ ‘(p) such that 
1/A(w) < €. Let us now construct a word ve Ml as follows: let 
a,, @ be generators of Jl,. Since p, # 0, a, occurs in w. Re- 


place one of the occurrences of a, in w by a,; denote the resulting 
word by v. Then we have 


A, (2) = d(x) -1 
A,(v) =A,(w) +1 
A (2) = A,(w), 3 < i) gn 


Since A(v) = A(w), we have 


A) A, (w) a ‘ ae 
3 - Mik ameiel ey | AC) 
Rte) AG) [| , is 

es ata) | 2G | 28? 

Cal ROU cart pan 

rayne ay a 


Hence ¢(v) lies in an €-neighborhood of p. 

Thus, the density of r(~) in a subset of &, is established. 
Since f is defined on r(~), we can ask whether or not f can be ex- 
tended to the closure of r(~) by the usual procedure of taking 
limits of sequences. If this extension can be performed, then the 
extended function thus obtained is continuous and uniquely deter- 
mined by 7. This is what we mean by the extension of f. 

We shall now show that the continuity of f (par abus de langage) 
implies the statement of Theorem 4. Let per(g) and let U(p, €) 
be an €-neighborhood of p. Then by Theorem 1 (Rosen, 1959b), it 
follows that p~*(U(p, €)) is a submonoid of {1 ,; in fact, it is the 
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quasi-ergodic submonoid which was denoted by E(p,, &). Pre- 
cisely the same reasoning implies that if ger(w) and V(qg, 5) is a 
S5-neighborhood of g, then ¥~*(V(q, 5) is a quasi-ergodic submon- 
oid of I,. If f is continuous, then the image of a connected set 
under f is mapped onto a connected set in #,; in particular, each 
set U(p, €) is mapped onto some V(q, 5), and the commutativity 
of the diagram on p. 200 implies precisely the statement of 
Theorem 4. 

On the basis of the foregoing, then, it is necessary to show that 
if P =(p,, Pa ees Pad pr = (p4’s Pass sees Pp, ) are close in Ro 
then f(p), f(p’) are close in #,. But we know from the previous 
calculation (Rosen, 1959a, p. 89) that 


- n 
a MiP i yy S;P i 
ACER ADA 1) Sp lea Sy P 
> + 8)P, Di (r+ 82, 
n n 


au. | | 492 Valen oe 
d+ ee, Yo + 827 
1 1 


Thus, we are concerned with the magnitudes of the quantities 


= TP i = Pi 
= (r, + 8,)p, x (1, + 8,)p; 


’ 


| = 8,P; 2X 8:P; 


E(r,+8,)p, U(r, +8,)p,|’ 2 


given that p, p’ are close, 
Now suppose that for each i, lp; —P; | <e& Then 


|Zrp,-Zrps| <re 
[2 sp,-2 sp7| <x se 
[E(r, + 8,)p,-X(r, + 8,)p7| <(20,+ 8) €. 


That is, specifying an order of closeness between p and p’ spe- 
cifies at the same time a degree of closeness between the respec- 
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tive numerators and denominators of the quantities (3) above. 
However, it is well known that the operation of taking quotients is 
a continuous operation (in any domain bounded away from zero); 
i.e. if |a—a’| and |o-’| are small, then | a/b — a/b’ is 
small. Applying this result to the quantities (3), we see that the 
continuity of f follows. This completes the proof. 

Let us now turn our attention to somewhat different matters. By 
way of motivation, we begin by drawing a simple corollary from 
theorems proved (Rosen, 1959a,b). Theorem 5 (1959a) asserts that 
every strictly ergodic submonoid of a finitely generated free monoid 
is free on a countably infinite set of generators. Theorem 3 (1959b) 
makes the same assertion for proper quasi-ergodic submonoids. For 
definiteness, let us suppose that EF is a strictly ergodic submonoid of 
tl,, E’ is a quasi-ergodic submonoid of {l,. Since E and E” are 
both countably generated, they are isomorphic (an isomorphism be- 
tween them being induced by any 1-1 correspondence between com- 
plete sets of generators of E and E’). Let f:E —» E’ be an iso- 
morphism. Then zt is impossible to extend f to a monomorphism of 
NM, into It, (i.e. it is impossible to find a monomorphism g :Il, — 
Ul, such that g(w) =f(w) for all we). For if such an extension 
were possible, it would follow from Theorem 6 (1959a) that g(E) is 
strictly ergodic in M,. By construction, g(£) = E” is not strictly 
ergodic in 2,; hence no such g exists. 

This result raises the question as to whether it is possible to 
generalize the above considerations as follows: if EF, E” are now 
taken to be guasi-ergodic submonoids of Ml,, Ul, respectively, is it 
possible to construct an isomorphism between E and E” which can- 
not be extended to a monomorphism of 9, into ,? We shall show 
that this is the case by constructing an example of such a non- 
extendable isomorphism. 

As usual, let {a,, Ayysees a} generate the free monoid | a and 
let io, b,} generate Il,. Let ECR, bea quasi-ergodic sub- 
monoid generated by the set W = iw, Wy yees 3, and let E’ C MN, be 
a quasi-ergodic submonoid generated by the set V= {0,5 Dey sini sia i 
In particular, we can find a word wv, eW and a word %%, eV such 
that A(w;.) > n(d, ). Let us now construct a 1-1 correspondence 
f between W and V such that f(w;,) = Cre By multiplicativity, f 
then extends to an isomorphism between E and E’. 

However, this isomorphism does not extend to a monomorphism 
f between M1, and l,. For if such an extension were possible, 
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then obviously each generator a, of I, must be mapped by f onto 
a word f(a,) of length greater than or equal to one; i.e. p [7(a;)] aus 
From this and the multiplicativity of f it follows that pf (w, I = 
A(u; ). On the other hand, we have by construction that 


wUf(w, 1 = w(v,.) <A(w,). 


This contradiction shows the impossibility of obtaining such an é 
Hence we have the result that, although a monomorphism f: 1,— 
ji, induces an isomorphism between a quasi-ergodic submonoid E 
of tt, and some quasi-ergodic submonoid E” of Il,, not every such 
isomorphism arises in this manner. 

Let us now investigate the implications of this result for the 
DNA-protein coding problem. We stated (1959a) that it was rea- 
sonable to suppose that any naturally occurring DNA-protein code 
was of an ergodic nature; i.e. was restricted to ergodic (or quasi- 
ergodic) submonoids of ‘l,, (the set of all polypeptide chains) and 
Ni, (the set of all DNA molecules), However, we explicitly re- 
stricted attention in loc. cit. to codes arising from restrictions of 
monomorphisms f:%,, > , to ergodic (or quasi-ergodic) sub- 
monoids of i,,, Il, respectively. For such codes, we showed that 
a finite set of conditions (namely, the stipulation of the numbers 
Pi» 7;, See especially Rosen, 1959a, pp. 91-92) was sufficient to 
completely derive the corresponding coding process, 

The result we have just obtained, however, makes it clear that 
our considerations in loc. cit. were too narrow; in particular, it is 
implied that there may exist admissible quasi-ergodic codes which 
cannot be derived from the restriction of monomorphisms /:Jl,, 
It, to suitable submonoids. Further, a moment’s reflection will re- 
veal that such ‘‘anomalous’’ codes cannot be completely deter- 
mined by any finite set of initial conditions, since quasi-ergodic 
Submonoids are countably generated. 

The possibility of the existence of ‘‘anomalous”’ DNA-protein 
codes in nature (and there seems to be no @ priort reason for ex- 
cluding them) considerably complicates the experimental biochem- 
ical picture. {n particular, the mathematical fact that these anoma- 
lous codes cannot be extended to the generators of Il,,, I, trans- 
lates into the biochemical statement that the study of amino acid 
and nucleotide compositions of polypeptides and DNA molecules 
related by such a code is completely irrelevant to the actual coding 
process. In fact, it is readily seen that the existence of anoma- 
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lous codes may contribute to the fact that all empirical correlations 
between amino acid composition of polypeptides and nucleotide 
composition of DNA attempted thus far have been unsuccessful 
(although the presence of ‘‘nonsense’’ DNA and the fact that many 
different quasi-ergodic coding processes of the orthodox type (i.e. 
derivable from monomorphisms f:%,, — Wl,) may be represented 
among the samples chosen for these correlations would also lead 
to such a result). 


This work has been aided by United States Public Health Serv- 
ice Grant RG-5181C2. 
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THE ROLE OF THE INDIVIDUAL IN SOCIAL DYNAMICS 
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The mathematical theories of social dynamics developed previously 
are deterministic and leave no room for the effects of single ‘‘excep- 
tional’? individuals. They all, however, lead to the existence of insta- 
bility points and threshold phenomena. It is shown that in the neighbor- 
hood of such instability points, exceptional individuals, who appear 
rarely, can appreciably advance or retard the moment at which the insta- 
bility is reached and at which a sudden change in the society occurs. 
Such individuals, however, do not either prevent or cause the instability, 
or the change, to occur. The indeterminacy introduced by ‘rare’? ine 
dividuals into the time course of social change is inversely proportional 
to the rate of social change. 


In our previous publications (for a summary, see Rashevsky, 
1947, 1959) we discussed different mathematical models of mech- 
anisms of social change. Those models describe the changes of 
behavior of a society with time. In all of them we find threshold 
phenomena; a behavior pattern of a society is stable up to a point, 
at which the smallest change of a parameter results in a relatively 
sudden change of the whole behavior pattern. The number of in- 
dividuals involved is one of the parameters. Thus the addition or 
the removal of a single individual at the critical point may pre- 
cipitate a finite change. 

Nevertheless all the above mentioned mechanisms are strictly 
deterministic and do not provide any special role in social change 
for a single individual. The individual whose addition to a group, 
or whose infinitesimal change in behavior precipitates a finite 
change in the behavioral configuration of the whole society, may 
be any individual. The situation is analogous to instable states 
in physics. Here also a displacement of a single molecule pre- 
cipitates a finite change of the physical configuration, but any 
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molecule can serve that purpose. There are no special individual 
molecules which produce physical changes in the system. Simi- 
larly, in the models mentioned above, there are no individuals who 
“make history.’’ Our mathematical theory of history does not bring 
in the single individual explicitly. 

Of course if the processes of social dynamics were determined 
completely and exclusively by single individuals, then no mathe- 
matical theory of social dynamics would be possible because of 
the unpredictability of the exact behavior of a single individual. A 
certain amount of determinism must be present in order that a math- 
ematical theory of social dynamics be possible. In our theories 
this determinism arises as a result of the interaction of a very 
large number of individuals, the behavior of each of whom is quasi- 
nondeterministic. Here we again have an analogy with physics. 

Yet observation indicates that sometimes single individuals do 
play an appreciable part in various social changes. Hence our 
models are incomplete in that they do not provide for such phe- 
nomena. This is not surprising, since our theories, like any theo- 
ries, are based on a large number of oversimplified assumptions. 
It is the purpose of this paper to show that a natural generalization 
of our previous considerations leads to a theory of the possible role 
of exceptional individuals in social dynamics. We assume fa- 
miliarity of the reader with our previous results (Rashevsky, 1959). 

Let us first consider the case of two active classes which ex- 
hibit mutually exclusive behavior, each of which attempts to in- 
fluence the third—passive—class to choose the corresponding 
behavior (Rashevsky, 1947, 1959, Chap. xiii). Let N, be the 
number of individuals in the society, X, and Y, the number of in- 
dividuals in the two active classes, @, and e, the corresponding 
coefficients of influence of the active individuals, and a, the co- 
efficient of influence of the passives. Then as we have seen, all 
passives will exhibit behavior R, of the class of X, individuals if 


C,—-a 
4 >—— Nj +—— Y, =a. (1) 
a + a, a+ a, 
If 
a a@,-a 
Y,>——-+_N, +2 * Xo, (2) 
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or, which is the same, 


Cy + a 


ae oe (3) 


0 
0 1 


then all passives exhibit behavior R, of the class of Yj) individ- 
uals. If X, is between the values of the right sides of equations 
(1) and (3), the passives may either all exhibit behavior R,, or all 
exhibit behavior F,. 

Let first inequality (3) be just satisfied; in other words, let X, 
be just slightly less than 8. All the passives then exhibit be- 
havior R,. To make the passives change their behavior to F,, X, 
must just exceed «. (As can be shown, we always have & > Ee) 
Hence X, must increase by the amount & — B. However, the sign 
of the inequalities (1) and (3) is also determined by the other pa- 
rameters entering into them, in particular by @, and ¢. For the 
sake of simplicity we shall consider only a. This does not in- 
troduce any limitations in principle. 


Putting 

aN +(€q -— 21) Vo - 21% sl (4) 

X 

0 

(¢, +a,)Y, -—a,N, + 4,% 

0 1 0 bey 6) 1250 se at (5) 

X 

) 

we may write inequalities (1) and (3) correspondingly as 

@>%13 4% < By: (6) 


It is readily seen that a, > 8, whenever X, + Y, < No, which is 
always true by definition. 

Now again, first let the second inequality (6) be just satisfied, 
so that all passives exhibit behavior F,. If then a, increases by 
an amount slightly exceeding 4, - B,, then the first inequality (6) 
will become satisfied, and the passives will all begin to exhibit 
behavior 2,. What can bring about such a change in @,? 

In our previous discussion we considered @, (as well as ¢, and 
a,) as fixed constants. The quantity @, is the coefficient of in- 
fluence of each individual of a class. However, even in a rather 
homogeneous class of people, no two individuals are alike in every 


respect. Therefore, the coefficient of influence @, actually rep- 
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resents an average over the X, individuals of a varying coefficient 
of influence a,, which is continuously distributed over X,. Let 
f(a) be the distribution function of a, so that f(a,)da, denotes 
the fraction of all individuals with an a, between a, and a, + dap. 


Then @, is given by: 
Z, -{ a,f(a,) da,. (7) 
0 


If f(a,) is defined in the range (—-, +), as are most distribu- 
tion functions, then, if X, were infinite, individuals with any a), 
as large as we wish, would be present. Actually, however, be- 
cause X, is not only finite but is also a positive integer, this is 
not so. The number of individuals with an a, > aj > 0 is given by: 


N (a) = X, | ; f(a) da,. (8) 
If 
‘. fe f(a,) da, << 1, (9) 


then there will not be in general, in any given generation, any in- 
dividuals with a, > a*. However, 


[,, 1(¢0) day = p (a8) (10) 


represents the probability of an individual having a, > a*. If gis a 
positive large number and if 


P (a>) = =< (11) 
4 

then about once in & generations an individual with a > af will 
appear. The larger ¢, the rarer the appearance of such an individ- 
ual, On the other hand, p (aj) decreases with increasing a*. Hence 
the larger a*, the rarer the appearance of an individual with ay > ap. 
While the average value @, of the coefficient of influence may 
remain relatively constant, yet when an individual with an ex- 
ceptionally large value a, > ay appears, the average @, may ap- 
preciably increase. If, as assumed, X, is sufficiently large, then 
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one individual with a high @, will change the average a, to 


q=-a,+—. (12) 


=) 
P< 


If 


ee B (13) 
senor Wee ce 
X 


then the appearance of such an individual will induce a change 
from behavior &, to behavior R,. The probability of such an event 
is given by p[X,(a«,-£,)], where p(x) is defined by equation 
(10). The larger 4, - 8,, the rarer such a possibility. For pre- 
scribed 4, and §,, which are given by equations (4) and (5), that 
is for prescribed characteristics of the society, the frequency of 
occurrence of such a phenomenon can be calculated if f(a,) is 
known. For large values of 4, —f,, such as are likely to occur, 
this frequency will be very small. 

Consider now the case where, at first, equation (3) is just satis- 
fied, so that the passives exhibit behavior R,. Let then X, and Y, 
change slowly until equation (1) becomes satisfied and the passives 
change to behavior #,. Such a slow change may, for example, be 
caused by the “thinning out”’ of Y,, as studied previously (Rashev- 
sky, 1947, chap. xiv; 1959, chap. xxiv). In such cases, inequality 
(1) begins to be satisfied at a certain time ¢*[cf. Rashevsky, 1947, 
p. 115, eq. (20); 1959, p. 209, eq. (20)], which is determined by the 
other parameters that characterize the society. In such changes 
we consider @, (and ¢,) as essentially constant, the effects now 
being due to changes of X, and Y,. Again, without loss of gen- 
erality, we shall consider the simpler case where only X, varies. 
At a time ¢*— Ad, the value of X, will be approximately, putting 
(dX ,/dt) 


t=t* = Vv; 


X,=4-vAt. (14) 


Since expression (1) is equivalent to the first inequality (6), at ¢* —- Az 
inequality (6) will not yet be satisfied. If, however, at that mo- 
ment an individual with a very large a* appears, then (6) may be- 
come satisfied. Therefore inequality (1) will also become satisfied, 
because even though X, has not yet reached the necessary value, 
@, has suddenly increased. 


212 N. RASHEVSKY 


With the assumption that only X, varies with time, Y, and other 
parameters remaining constant, we see that & remains constant, as 
long as an exceptional individual does not appear. 

At ¢* - Ad the value a,(¢* - AZ) of a, is given, because of equa- 


tions (4) and (14), by 


a,N, +(é, -4,)Y¥, 2 


ee = ; 15) 
a, (¢*-Ad) ee . ( 
At t* we have X, = a; therefore, because of (14), 
a (¢* — Az) > a, (¢*). (16) 


In order for the change of behavior to occur at ¢* — AZ instead of at 
¢*, an exceptional individual must appear at ¢* —~ At, who will in- 
crease @, according to (12) to such an extent that 


z > a,(¢*- Az) - a, (¢*) (17) 
or, because of (14): 
a, >(a -vAt) [a,(¢*- Az) - a ,(¢*)]. (18) 
Because of equations (4), (14), and (15) this may be written: 


. aN, +(e, -a,)¥, 


vAt. (19) 
a 


The probability of occurrence of inequality (19) decreases with in- 
creasing Az. It is again given by equation (10), into which we 
introduce (19). 

Instead of the appearance of an individual in the X, class with 
a very large a), we can consider in a similar manner the appear- 
ance in the Y, class of an individual with a very large c,. This 
will result in a delay by A¢ of the transition to behavior R . 

If, as we have done elsewhere, we consider both the variation of 
X, and Y, (Rashevsky, 1947, 1959), we obtain more complicated 
expressions but the principle remains the same. 

We find that the appearance of exceptional individuals intro- 
duces an uncertainty A¢ in the time of the sudden change of be- 
havior. It does not, however, prevent the occurrence of this change. 

We notice that the required value of a) increases also with »v, 
the rate of change of the social configuration. The rapidly oc- 
curring, short range sociological or historical processes are there- 
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fore less likely to be affected by rare individuals than are long 
range processes, 

Moreover, for the same value of a, that is for the same probabil- 
ity, we have from (19), denoting by C the fraction on the right side: 

C 
At=—. (20) 
v 
Thus for a given probability of the uncertainty, the latter is in- 
versely proportional to the rate of the social process. 

As we have seen elsewhere (Rashevsky, 1960), the equations of 
social interaction which lead to equations (1) and (2) are first- 
order approximations to a more elaborate theory which involves 
neurobiophysical factors. The latter expresses the change w in 
the tendency ¢ of every individual toward Ff, or F,, in terms of 
the number X of passives which at any given moment exhibit be- 
havior # , and the number Y of passives which exhibit behavior F,: 

dy 
te ee ayo. 8o, (21) 
where A and a are constants. 

The quantity X — Y considered as a function of yw is represented 
by a sigmoid curve. If the distribution function N(¢) of ¢ in the 
population is symmetric with respect to ¢=0, then the sigmoid 
curve passes through y = 0 and is symmetric with respect to the 
origin. Equation (21) has in general three solutions, one unstable 
and two stable. They are given by the intersections of the sigmoid 
curve with the straight linea J + @,X, -@,Y,. For @X) - ¢)¥ = 
0, the straight line passes through the origin, and there are two 
stable solutions. The effect of the term @,X,-c¢,Y, is to shift 
the sigmoid curve with respect to the straight line along the y- 
axis by the amount (@,X, — €,Y,)/a. The amount of the shift nec- 
essary to bring the society from one behavioral configuration, F,, 
to another, &,, is given approximately (Rashevsky, 1960) by 


2 ABN, - a)® 
ages ee (ae stella (22) 
Le Sea hor tee 


Hence if originally the society exhibited behavior F,, but if the 
inequality a - ni a ae ee Dh 
GX yal 6 > Oo: OF. 4 > mo a 
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becomes satisfied, the behavior pattern shifts to R,. Since ini- 
tially we had @X, = ¢)Y,, or @ =¢)Y,/X,, the total change of 
@, must be equal to 6,/X,, which according to equation (12) im- 
plies the appearance of an individual with an a > 6,. The prob- 
ability of such an event is given by p(d,), where p(z) is defined 
by equation (10) and 6, is given by equation (22). All the previous 
considerations remain unchanged except for the substitution of 
6, fora, - B,.- 

When X, = ¥, = 0 in the more general theory, we have a purely 
imitative behavior. In this case there is no effect of the individual 
at all. It might be thought that individuals with very large values 
of ¢, which will appear once in a while as determined by N(¢), 
would produce that effect. It must be remembered, however, that 
an individual with ¢ =+e merely exhibits always behavior RE, 
while one with ¢ =—« exhibits always behavior R,. If No, and 
therefore X and Y are very large, then the appearance of one or 
even a few such individuals affect X — Y quite insignificantly. 

The situation would change if we assume, as has been done in 
a previous paper (Rashevsky and Householder, 1941; Rashevsky, 
1947), that the higher the absolute value of ¢ for a given individ- 
ual, the greater his coefficient of influence in the given direction. 
The above-mentioned paper does not however provide a sufficient 
framework for such a theory because the coefficient of influence is 
assumed there to remain always finite. To obtain the desired ef- 
fect, the coefficient of influence should increase to infinity for 
@=t+o, Such cases present interesting features and should be 
made the object of a separate study. 

It may be remarked that the assumption that a person with a 
Strong tendency toward a given behavior also has a large coef- 
ficient of influence does not seem to be corroborated by facts. 
An individual may be obsessed with preaching some unusual doc- 
trine yet may remain without any influence on others. 


This work was aided by a grant from the Dr. Wallace C. and Clara 
A. Abbott Memorial Fund of the University of Chicago. 
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